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Abstract 

Wavelets are a useful basis for constructing solutions of the inte- 
gral and differential equations of scattering theory. Wavelet bases effi- 
ciently represent functions with smooth structures on different scales, 
and the matrix representation of operators in a wavelet basis are well- 
approximated by sparse matrices. The basis functions are related to 
solutions of a linear renormalization group equation, and the basis 
functions have structure on all scales. Numerical methods based on 
this renormalization group equation are discussed. These methods 
lead to accurate and efficient numerical approximations to the scat- 
tering equations. These notes provide a detailed introduction to the 
subject that focuses on numerical methods. We plan to provide peri- 
odic updates to these notes. 
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1 Introduction 

Wavelets are versatile functions with a wide range of applications includ- 
ing time-frequency analysis, data compression, and numerical analysis. The 
objective of these notes is to provide an introduction to the properties of 
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wavelets which are useful for solving integral and differential equations by 
using the wavelets to represent the solution of the equations. 

While there are many types of wavelets, we concentrate primarily on 
orthogonal wavelets of compact support, with particular emphasis on the 
wavelets introduced by Daubechies. The Daubechies wavelets have the ad- 
ditional property that finite linear combinations of the Daubechies wavelets 
provide local pointwise representations of low-degree polynomials. We also 
have a short discussion of continuous wavelets in the Appendix I and spline 
wavelets in Appendix II. 

There notes are not intended to provide a complete discussion of the sub- 
ject which can be found in the references given at the end of this section. 
Rather, we concentrate on the specific properties which are useful for nu- 
merical solutions of integral and differential equations. Our approach is to 
develop the wavelets as orthonormal basis functions rather than in terms of 
low- and high-pass filters, which is more common for time-frequency analysis 
applications. 

The Daubechies wavelets have some properties that make them natural 
candidates for basis functions to represent solutions of integral equations. 
Like splines, they are functions of compact support that can locally point- 
wise represent low degree polynomials. Unlike splines, they are orthonormal. 
More significantly, only a relatively small number of wavelets are needed to 
represent smooth functions. 

One of the interesting features of wavelets is that they can be generated 
from a single scaling function, which is the solution of a liner renormalization- 
group equation, by combinations of translations and scaling. This equation, 
called the scaling equation, expresses the scaling function on one scale as a 
finite linear combination of discrete translations of the same function on a 
smaller scale. The resulting scaling functions and wavelets have a fractal-like 
structure. This means that they have structure on all scales. This requires a 
different approach to the numerical analysis, which is provided by the scaling 
equation. These notes make extensive use of the scaling function. 

Some of the references that we have found useful are: 
[1] I. Daubechies, Orthonormal bases of compactly supported wavelets, Comm. 
Pure Appl. Math. 41(1988)909. 

[2] G. Strang, "Wavelets and Dilation Equations: A Brief Introduction," 
SIAM Review, 31:4, pp. 614-627, (Dec 1989). 

[3] I. Daubechies, Ten Lectures on Wavelets, SIAM, Philadelphia, 1992. 

[4] C. K. Chui Wavelets - A tutorial in Theory and Applications, Academic 
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Press, 1992. 

[5] W.-C. Shann, "Quadrature rules needed in Galerkin-wavelets methods", 
Proceedings for the 1993 annual meeting of Chinese Mathematics Associa- 
tion, Chiao-Tung Univ, (Dec 1993). 

[6] W.-C. Shann and J.-C. Yan, "Quadratures involving polynomials and 
Daubechies' wavelets", Technical Report 9301, Department of Mathematics, 
National Central University, (1993). 

[7] G. Kaiser, A Friendly Guide to Wavelets, Birkhauser 1994. 
[8] W. Sweldens and R. Piessens, "Quadrature Formulae and Asymptotic 
Error Expansions for wavelet approximations of smooth functions", SIAM J. 
Numer. Anal., 31, pp. 1240-1264, (1994). 

[9] H. L. Resnikoff and R. O. Wells, Wavelet Analysis, The Scalable Structure 
of Information, Springer Verlag, NY. 

[10] O. Bratelli and P. Jorgensen, Wavelets through a Looking Glass, Birkhauser, 



In addition, some of the material in these notes is in our paper 

[11] B. M. Kessler, G. L. Payne, and W. N. Polyzou, Scattering Calculations 
With Wavelets, Few Body Systems, 33,1-26(2003). 

2 Haar Scaling Functions and Wavelets 

Scaling functions play a central role in the construction of orthonormal bases 
of compactly supported wavelets. The scaling functions and wavelets are 
distinct bases related by an orthogonal transformation called the wavelet 
transform. 

The concept of scaling functions is most easily understood using Haar 
wavelets (these are made out of simple box functions). The Haar functions 
are the simplest compactly supported scaling functions and wavelets. 

The Haar scaling function is defined by 



2002. 




x < 

1 < x < 1 
x > 1 



(1) 



It satisfies the normalization conditions: 

POO pi 




(0,0):= / (jf {x)4>{x)dx = I 4>{x)dx = 1 



(2) 
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The operations of discrete translation and dilatation are used extensively 
in the study of compactly supported wavelets. The unit translation oper- 
ator T is defined by 

(T X )(x)= X (x-l). (3) 

This operator translates the function x( x ) to the right by one unit. The unit 
translation operator has the property: 

/oo 
r(x-l)x(x-l)dx= (4) 
-oo 



oo 

oo 



ii*(y)x(y)dy = (^,x) (5) 

> 

where y — x — 1. This means that the unit translation operator preserves 
the scalar product: 

(T^,Tx) = ^,x)- (6) 
If A is a linear operator its adjoint is defined by the relation 

ty,At X ) = (Ail;,x). (7) 

It follows that 

/oo 
r(x-l) X (x)dx. (8) 
-oo 

Changing variables to y = x — 1 gives 

/oo 
r(y)x(y + i)dy (9) 
-oo 



or 



{T\){x) = X {x + l) (10) 
which is a left shift by one unit. Since 

(^,x) = (T^T X ) = (^TtT X ) (11) 

it follows that T" 1 " = T _1 . An operator whose adjoint is its inverse is called 
unitary. Unitary operators preserve inner products. 

It follows from the definition of the Haar scaling function, <f>(x), that 

/oo 
(j)*(x)(j)(x-n + m)dx = 
-oo 



1 

<f>(x — n + m)dx = 5 nm (12) 

o 

This means the functions 

4> n {x):={T n <f>){x) = 4>{x-n) (13) 

are orthonormal. There are an infinite number of these functions for integers 
n satisfying — oo < n < oo. 

The integer translates of the scaling function span a space, Vo, which is a 
subspace of the space of square integrable functions. The elements of Vo are 
functions of the form 

oo oo oo 

/(*)= E fnM%)= E /n(T»(x)= E fn<l>(x-n), (14) 
n=— oo ra=— oo ra=— oo 

where the square integrability requires that the coefficients satisfy 

oo 

E \fn\ 2 < OO. (15) 
n=— oo 

For the Haar scaling function Vo is the space of square integrable functions 
that are piecewise constant on each unit-width interval. Note that while 
there are an infinite number of functions in Vo, it is a small subspace of the 
space of square integrable functions. 

In addition to translations T, the linear operator D, corresponding to 
discrete scale transformations, is defined by: 

(Dx)(x) = -j=x(x/2). (16) 

When this is applied to the Haar scaling function it gives 

{0 x < 
73 0<x<2 . (17) 
x > 2 

This function has the same box structure as the original Haar scaling 
function, except it is twice as wide as the original scaling function and shorter 
by a factor of Note that the normalization ensures 

/■oo 1 

(D^,D X ) = J -r(x/2) X (x/2)dx (18) 
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r(y)x(y)dy = (Ax) (19) 



where the variable in the integrand has been changed to y — x/2. 
The adjoint of D is determined by the definition 

/OO 1 
-^=r(x/2)x(x)dx. (20) 



Setting y = x/2 gives 



which gives 



/OO 
r(y)V2 X (2y)dy (21) 
-oo 



(D^ X )(x) = V2 X (2x). (22) 

This shows that = D^ 1 or .D is also unitary. 

Define the functions constructed by n translations followed by m scale 
transformations 

cf> mn {x) = (D m T n (j))(x) = {D m (f) n ){x) (23) 

= 2~ m/2 <p{2~ m x -n) = 2- m/2 <p{2- m (x - 2 m n)). (24) 
It follows that for a fixed scale m 

(<f>mn, <p mk ) = (D m <f> n , D m <p k ) = (<f> n , D m ~ m (p k ) = (0 n , <t> k ) = 5 nk . (25) 

This shows that the functions <f) mn {x) for any fixed scale m are orthonormal. 

We define the subspace V m of the square integrable functions to be those 
functions of the form: 

OO OO 

fix) = Mmnix) = UD m T n 4>){x) (26) 

n=— oo ra=— oo 

where the square integrability requires that the coefficients satisfy 

oo 

E < °°- ( 27 ) 

n=— oo 

These elements of V m are square summable functions that are piecewise con- 
stant on intervals of width 2 m . The spaces V m and Vo are related by m scale 
transformations D m Vo = V m . 



6 



In general the scaling function 4>(x) is defined as the solution of a scaling 
equation subject to a normalization condition. The scaling equation relates 
the scaled scaling function, (D<f))(x), to translates of the original scaling 
function. The general form of the scaling equation is 

(D<f ) )(x) = J2hiT l <J>(x) (28) 
i 

where hi are fixed constants, and the sum may be finite or infinite. This 
equation can be expressed as 



^(f) = EM*-*) ( 29 ) 



1 , ,x. 

i 

which is sometimes written as 

(f>(x) = V2^2hi(j)(2x-l) =^2ci<l)(2x-l) (30) 
i i 

where q = V2hi. Equation (|28j) is the most important equation in these 
notes. 

In general the scaling equation cannot be solved analytically. In the 
special case of the Haar scaling function the solution is obtained by observing 
that the scaled box is stretched over two adjacent boxes with a suitable 
reduction in height. It follows that: 

D<f>(x) = -j=<t>(x/2) = ±=cj>{x) + ^fcj>{x) 

= * 0(x) + i=#r-l). (31) 

Here ho — h\ — l/y/2. These coefficients are special to the Haar scaling 
function. The best way to think about the scaling function <f>(x) is to note 
that the scaling function <p{x) is the solution of the scaling equation up to 
normalization. The normalization is fixed by 

4>(x)dx = 1. 

An additional relation involving the translation T and dilatation operator 
D is useful for future computations. First note that 
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DT^x) = D^(x - 1) = - 1) = -jfC-^) = T 2 D^j(x), (32) 

which leads to the operator relation 

DT = T 2 D. (33) 

It follows from this equation that 

D(f) n (x) = DT n <f){x) = T 2n D</){x) = T 2n {h <f){x) + h{r<j>{x)). (34) 

This shows that all of the basis elements in V\ can be expressed in terms of 
basis elements in V . For the case of the Haar scaling function this is obvious, 
but the argument above is more general. 
Specifically if tp(x) G Vi then 

oo oo 
= ^Inix) = ^ d nD<p n (x) (35) 

n=— oo n=— oo 

oo oo 

= ^2 [d n h (i)2n(x) + dnhfon+iix)] = ^2e n (j) n (x) (36) 

n=— oo — oo 

where 

&2n — d n ho &2n+l = d n h\. (37) 

It is easy to show that 

oo oo 
n=— oo n=— oo 

What we have shown, as a consequence of the scaling equation, is the 
inclusion property 

V D Vi. (39) 
Similarly, using the same method, it is possible to show the chain of inclusions 

• • • V_ fc D V_ fe+1 D • • • D V D • • • V k D V k+1 ■ ■ ■ (40) 

These properties hold for the solution of any scaling equation. In the Haar 
example the spaces V rn are spaces of piecewise constant, square integrable 
functions that are constant on intervals of the real line of width 2 m . 
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The subspaces V m are used as approximation spaces in applications. To 
understand how they are used as approximation spaces note that asm^ — oo 
the approximation to f(x) given by 

oo 

fm(x) = ^ fmn<Pmn(x) (41) 
n=— oo 

with 

/oo 
<p mn {x)f{x)dx (42) 
■oo 

is bounded by the upper and lower Riemann sums for steps of width 2 _m . 
This is because, up to a scale factor, the coefficients f mn are just average 
values of the function on the appropriate sub-interval (to deal with the infinite 
interval it is best to first consider functions that vanish outside of finite 
intervals and take limits). Since the upper and lower Riemann sums converge 
to the same integral (when the function is integrable) it follows that 



J — c 



|/ m (aO-/(a;)|tfa;<e (43) 



for sufficiently large —to. A similar argument can be extended to get L 2 
convergence . 

Similarly, as to — > +00, the width of 4> mn (x) grows like 2 m while the height 
falls off like 2" m / 2 . A gain, if the function vanishes outside of a bounded 
interval then for sufficiently large to there is only one (or two) <p m n{x) that 
are non-vanishing where the function is non-vanishing. In the case that only 
one (p mn = (fimno overlaps the support of f(x) 

/oo 
f(x)dx. (44) 
-00 

The integral of the square of this function ~ 2~ m — > as m — > 00. 
Note that 

/oo poo 
f m (x)dx^ / f(x)dx (45) 
-00 J —00 

as m — > 00. This shows that the limit of the integral of f m {x) as m — > 00 is 
finite in L 1 but in L 2 . 
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It is useful to express some of these results in a more useful form. Define 
the projection operators 

oo 

(x) (46) 

n=— oo 

where 



fmn = / (j)* mn (x)f(x)dx. (47) 

J — oo 

The above conditions can be stated in terms of these projectors: 

lim P m = I (48) 



m— ►— oo 



lim P m = 0. (49) 

m— >+oo 

These results mean that the approximation space V m approaches the space 
of square integrable functions as m — > — oo. We have shown that (|48|) and 
f!49|) are valid for the Haar scaling function, but they are also valid for a large 
class of scaling functions, 

We are now ready to construct wavelets. First recall the condition 

V D Vi. (50) 

Let Wi be the subspace of vectors in the space Vo that are orthogonal to the 
vectors in V\. We can write 

Vo = VieWi. (5i) 

This notation means that any vector in Vo can be expressed as a sum of two 
vectors - one that is in V\ and one that is orthogonal to every vector in Vi. 

Note that the scaling equation implies that every vector in V\ can be 
expressed as a linear combination of vectors in Vo using 

D<p n (x) = h Q (p2n(x) + /li0 2 „+l(x). (52) 

Clearly the functions that are orthogonal to these in Vi on the same interval 
can be expressed in terms of the difference functions 

i/ji n (x) := Di/; n (x) = hx4> 2n {x) - ho4> 2 n+i(x) = -^=(4> 2n (x) - 4> 2n +i{x)). (53) 
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Direct computation shows that the ip\ n {x) are elements of Vo that satisfy 

(Dipi n , D<pi) = 0. (54) 

and 

(ipln^lk) = Kk- (55) 

Thus we conclude that Wi is the space of square integrable functions of 
the form 

oo 

fix) = Mm(x) (56) 

n=— oo 

with 

oo 

/(*) = E i/«i 2 - ( 57 ) 

n=— oo 

Similarly, we can decompose Vi = V;+i © Wi+i for each value of I. For 
the special case of W we define the Haar mother wavelet as 

iP(x) := D-^h^x) - h T(f)(x)) = (58) 

/iiv / 20(2t) - h V2(/>(2(t - 1)) = (0(2t) - 4>{2{t - 1))) (59) 

which is manifestly orthogonal to the scaling function. Translates of the 
mother wavelet define a basis for W 

t/j n (x) = T>(x) = T n D' l (h 1( f>(x) - hoT(f>(x)) = (60) 

D^ihfonix) - ho<p2n+l{x)). (61) 

If we decompose V m we have: 

V_ m = W_ m+1 © V_ m+ i 

= W_ m+ i © W- m+ 2 © v_ m+2 

© yv-m+2 © • • • ©Wj©V t . (62) 

Note that unlike the V m spaces, the W m spaces are all mutually orthogonal, 
since if m > n — > W m C V n which is orthogonal to W n by definition. 
If f(x) is any square integrable function the conditions 
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lim P m = I (63) 

TO 
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lim P m = (64) 

m— >+oo 

mean that for sufficiently large m and any / that f(x) can be well approxi- 
mated by a function in 

W_ m+ i © W_ m+2 © ■ ■ • © W,. (65) 

This means that the function can be approximated by a linear combination 
of basis functions (wavelets) from each of the spaces W r . 

A multiresolution analysis is a set of subspaces V m and W m satis- 
fying (|SHj), and The condition (j6Hj) allows one to interpret the 
space V m , for sufficiently large —m, as an approximation space for numerical 
applications. 

Basis functions for W m are given by 

i, mn {x) = D m T n ^{x) = D m -\h x <t> 2n {x) - Man+i(aO). (66) 

That these are a basis with the required properties is easily shown by showing 
that these functions are orthogonal to V m and can be used to recover the 
remaining vectors in V m -i- 

The functions ip n i(x), are called Haar wavelets. They satisfy the orthonor- 
mality conditions: 

(i>rd,i>n'l f ) = Srm'SiV (67) 

where the 5 nn i follows from the orthogonality of the spaces W n and W n > for 
n 7^ n'. 

The Su> follows from the unitarity of D and 

(^,T»=5 n0 . (68) 

The important steps discussed above generalize to the case of a general 
scaling equation of the form: 

D(j){x) = ^h{T l 4){x). (69) 

This equation is solved to find the scaling function 4>{x). This, along with 
translations and dilatations is used to construct the spaces Vi. The scaling 
equation ensures the existence of spaces W m , satisfying V m+ \ = W m © V m 
that can be used to build discrete orthonormal bases. The mother wavelet 
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function is expressed in terms of the scaling function and the coefficients hi 
as 

^x) = D- l Y,9iT l <P(x) (70) 
i 

where we will see later that 

9i = {-) k h k ~i (71) 

where k is any odd integer. In general the coefficients hi must satisfy con- 
straints for the solution to the scaling equation to exist. General wavelets 
can be expressed in terms of the mother wavelet using ([66)1 . In the next 
section the coefficients gi will be expressed in terms of the scaling function. 

3 Scaling Functions - General Considerations 

This section extends the treatment of scaling equation to a more general class 
of scaling functions than the Haar functions. In general, a scaling function 
satisfies the following three conditions. First, the scaling function is the 
solution of the scaling equation 

D<f>(x) = ^2hiT l (f>(x) (72) 
i 

where hi are numerical coefficients that define the scaling equation. Second, 
in addition to satisfying the scaling equations, integer translates of the 
scaling functions are required to be orthonormal 

(<f> n , m ) = (r>, r» = (</>, T m -» = s mn . (73) 

Third, the initial scale is fixed by the normalization condition 

/«,)*- 1. (74) 

It might seem like the normalization conditions in (|T3*j) and (J71j) are not 
compatible. To see that this is not true note that condition (JT3*|) is invariant 
under unitary changes of scale of the form 



D s x(x) := 
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while condition (J73J) is not. It follows that condition (J73J) can be interpreted 
as setting a starting scale, s. The condition ()73|) is preserved independent of 
the starting scale. 

We now investigate the consequences of these three conditions. Using the 
definitions of the operators D and T the scaling equation becomes: 

^(f) = £M(*-l). (75) 

As shown in section 1, it can be put in the useful form 

<f>(x) = ^2V2ht<p(2x-l). (76) 
i 

In general the sums may be from — oo — ► oo. Finite sums are treated by 
assuming that only a finite number of the h^s are non zero. All of the 
compactly supported scaling functions are solutions of scaling equations with 
a finite number of non-zero coefficients. 

If the scaling equation has a solution, it is unique up to an overall nor- 
malization factor. To see this take the Fourier transform of both sides of 
equation (J77)j) to get 

1 poo 1 poo 

4>(k) = -= / e~ lkx <p(x)dx = V V2ht-= / e~ ikx <P(2x - l)dx. (77) 

V27T J-cg t V27T J-oo 

Changing variables x — > 2x — I on the right-hand side gives 

-I /*oo 1 1 /*oo 

-= / e- ikx (j)(x)dx = V -^1-= / e- i(fc/2){a;+/) 0(a;)rfx (78) 

V27T J-oo ^ V2 V27T J-oo 

or 

= Q) * (|) (79) 

where 

~ h (k) = E 7= e " ifc • ( 8 °) 
This form of the scaling equation can be iterated n times to get: 

^ ' m=l ^ ' 
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This equation holds for any n provided the Fourier transforms exist. For 
a finite n, an approximation can be made by a finite number of iterations of 
the form 



Mk) = <t>n-i ^- ) h ^- j (82) 

for any starting function <f>o(k). In the limit of large n the function (f> n (k) 
should converge to a solution to the scaling equation. The result of formally 
taking this limit is 

to) = 2£.te*)iw?) 

1=1 

oo , 

= mHH- l ). (83) 
i=i 

If the limit exists as n — > oo, and the scaling function is continuous in a 
neighborhood of zero, then the solution of the scaling equation is uniquely 
determined by the scaling coefficients hi up to the overall normalization <f>o(Q). 
The condition 0o(O) = l/y/2n is equivalent to the standard normalization 
condition 

roc 

4>(x)dx = 1. 



j — c 



The resulting solution of the scaling equation is independent of the choice 
of starting function provided it is normalized so 0(0) = l/v / 27r. Once the 
normalization is fixed, the limit only depends on the coefficients hi. 

Thus, if the infinite product converges, then we have an expression for 
the scaling function, up to normalization, which is fixed by assigning a value 
to 0(0). To show how this works we compute this limit for the Haar scaling 
equation. 

For the Haar scaling equation the expression for the scaling function is 



tlHa+e-^) 



= lim -L4(l + e~ tk/2 )(l + e~ ik ' A ) ■■■{! + e"^ 2 ' 
expanding this out in powers of e~ tk ^ L gives 

2 l -l 

1 limi V( e - ifc / 2; r 



V2n *-°o 2 l ^ n 

v m=0 



15 



1 1 1 - e~ ik 

lim 



A^f i^oo 2 l l- e~ ik / 2 ' 



1 „-ifc/2 sin ( fc / 2 ) 



v 7 ^ (k/2) K ! 

A direct calculation of the Fourier transform of the Haar scaling function 
gives 

1 Z" 1 



(f)(k) = —= / e~ lkx (j)(x)dx = -= I e~ lkx dx 
y2ix J-oo V27T 

_ 1 .-ifc/2 sin ( fc / 2 ) 



^ {k/2) (85) 

which agrees with (|84|) . 

The above analysis shows that the solution of the scaling equation de- 
pends on the choice of scaling coefficients hi. The scaling coefficients hi are 
not arbitrary. First note that setting k = in (J83|) gives 



H~h(0). (86) 



1=0 



Now using (|50|l gives 



MO) = i = E 7= ( 8? ) 



or 



^2hi = V2. (88) 
z 

This condition is satisfied by the Haar wavelets. This is a necessary condition 
on the scaling coefficients in order to have a solution to the scaling equation. 

Another condition which constrains the scaling coefficients is the orthog- 
onality of the unit translates, (<f> n ,<f> m ) — bnm- This requires, using f7fij) . 

/oo 
(f)(2x -2n- l)(f)(2x - 2m - k)dx 

/oo 
<p(2x)<p(2x - 2(m -n)-(k- l))dx 
■°° 

/oo 
4>{x)<j){x — 2(m — n) — (k — l))dx 
■oo 



16 



— hlhi-2(m-n) — <5mn (89) 
I 

or equivalently 

^hi^mhi = <W (90) 
i 

This is trivially satisfied for the Haar wavelets. Here and in all that follows 
we restrict our considerations to the case that the scaling coefficients and 
scaling functions are real. 

The orthogonality condition also requires that the number of non-scaling 
coefficients must be even. To see this assume by contradiction that there are 
2N + 1 non-zero scaling coefficients, ho ■ ■ ■ h 2 N- Then setting m = —N in 
(HHJ) gives 

22 h i+2Nhi = h 2N h = 5 N0 = 0. (91) 
i 

which means that either h = or h 2 N = 0, which contradicts the assumption 
that there are 2N + 1 non-zero scaling coefficients. This shows that if the 
number of non-zero scaling coefficients are finite, then there must be an even 
number, 2K, with / = • • • 2K — 1. 

Note that if there are only two non-vanishing scaling coefficients, ho and 
hi, then the conditions (|88|) and (J9U|) have a unique solution, which is the 
Haar scaling coefficients. In this case these equations become 

ho + h l = V2 (92) 

h h + hihi = 1. (93) 

These equations have the unique solution h = hi = 1/ \[2. 

Conditions ()88|) and ()90|) are important constraints on the scaling coeffi- 
cients. 

For scaling equations with more than two non-zero scaling coefficients, 
additional conditions are needed to determine the scaling coefficients. 

The number of non-zero scaling coefficients determines the support of 
the scaling function. The important property is that scaling functions that 
are solutions of a scaling equation with a finite number of non-zero scaling 
coefficients have compact support. The support is determined by the number 
of non-zero scaling coefficients. 

To determine the support of the scaling function, consider a scaling equa- 
tion with N = 2K + 1 non-zero scaling coefficients. The scaling function is 
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given by 

m = S f e ' fa n 

* 00 m=l 

iV-1 iV-l / m , \ m 

=J»E-E II Ewn. 

ni=0 n m =0 \fe=l v / fc=l 

This defines the scaling function as a distribution. This is not a useful 
representation for computation, however it indicates that if a scaling function 
has N non-zero coefficients hi then the scaling function has support on 

[0,{N-l)(\ + \ + \-)] = [0,N-l] 

where N is the number of non-zero scaling coefficients. 

While the support condition depends only on the number of non-zero coef- 
ficients, there are many scaling functions with N non-zero scaling coefficients. 
Except for the constraints dictated by the scaling equation, orthonormality, 
and normalization, there is considerable freedom in choosing the coefficients 
hi. 

The scaling coefficients also determine the mother wavelet function. In 
the general case the spaces V m are the spaces of square integrable functions 
spanned by the orthogonal basis functions <f> mn {x) := D m T n (p(x) for integer 
n satisfying — 00 < n < 00. As in the Haar case, the scaling equation implies 
that V m 3 V m +k for k > 0. Wavelet spaces are defined by 

W m : Vmr-l = V m © W m . 

It they also satisfy and they define a multiresolution analysis. 
The mother wavelet function lives in the space Wo which means that it 
has an expansion in V_i: 

if>(x) = y/29n<K2x - n) = ^^^r^x). (94) 

n n 

This equation can be expressed in a form similar to the scaling equation: 

Di/;(x) = J29nT n (f>(x). (95) 
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The mother wavelet and all of its integer translates should be orthogonal 
to the scaling function, which is in Vo- In terms of the coefficients this 
requirements is: 

(V> m , <fi) = 2j hig n ((p n+2m , (pi) 

n.l 



^ hg n 8, 



n+2m,l 



n.l 



= 2j K+2m9n = (96) 
n 

for all m. 

Orthonormality of the translated mother function requires 

(lpm,^n) = y^ J 9l9k((fil+2m,(fik+2n) 

l,k 

2_j9k+2(n-m)9k — 8 mn (97) 
k 

or equivalently 

(VVn, = 2J 9k+2m9k = <W (98) 

The choice 5^ : = (—{) h hi-k where I is any odd integer it satisfies (|93|) and 
(HU): 

9k+2(n-m)9k = ^^( — l) fe+2 ^ n ^ hi_ k _2(n~m) ( — l)*^J-Jfc 

= hk>+2(n-m)hk> = 5 mn (99) 

fe' 

where we have let fc' = / — /c in the last term. It also follows that 

2_. hn+2m9n = ^n+2m(~l)"^l-n 

n n 

= hl- n '(— I) 1 n 2m h n i + 2m = ( — )' hi^ n >{ — 1)™ h n i + 2m 

n' n' 

= (-l) l J29n'hn>+2m. (100) 
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Since I is odd, the sum is equal to its negative which shows that it vanishes. 
The choice of I is arbitrary - changing I shifts the origin of the mother by an 
even number of steps. Since the mother is orthogonal to all integer translates 
of the scaling function, it does not matter where the origin is placed. 
This shows that the coefficients hi, satisfying 

J2 h i = ^- ( 101 ) 
i 

hi-2mhi = S m0 (102) 

i 

with g k defined by 

g k := (-l) k hi- k I odd (103) 

give a multi-resolution analysis, scaling function, and a mother function. 
The Daubechies order- K wavelets are defined by the conditions 

= -M, (104) 

These equations ensure that polynomials of degree < K — 1 can be lo- 
cally represented by finite linear combinations of scaling functions on a fixed 
scale. This is a useful property for numerical approximations. The order 
i^-Daubechies scaling function has 2K scaling coefficients, with K = 1 cor- 
responding to the Haar wavelets, and each additional value of K adds one 
more orthogonality condition. 

The scaling equation (J95)) and the moment conditions (|104|) for the mother 
wavelet function gives the additional equations necessary to find the Daubechies 
scaling coefficients, hi. 



= (x n ,^) = (Dx n ,D^) 
J dxx n 2- n - 1,2S ^g. 



mrK x - m) 

m 



This gives 

/ dx(x + m) n g m (p(x) = 0. 

m J 

For n = this gives (using the n = equation) 

X> m = o,^x;(-iHi,_ ra = o, 
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Table 1: Scaling Coefficients 



hi 


K=l 


K=2 


K=3 


h 


1/V2 


(1 + V3)/aV2 




v/T0+ \/5 + 2y 




h 


1/V2 


(3 + VS)/AV2 


(5 + ^ 


v / T(D + 3v / 5 + 2 


v/10)/16-\/2 


h 2 





(3 - y/3)/4y/2 


(10- 


2^ + 2^5 + 


2 v / T0)/16v / 2 


h 3 





(1 - Vs)/aV2 


(10- 


2^-2^5 + 


2^)716^2 


/i 4 








(5 + ^ 


/TO -3^5 + 2 


v / 10)/16v / 2 


h 5 










/TO - V5 + 2V 


To)/i6-\/2 



for n = 1 this gives 

mg m = 0,^J2 m{-l) m hi. m = 0, 

m 

for n = 2 • • • this gives 

^ m 2^ = 0, ^ ^ m 2 (-l) m /^ m = 0, 



^ m ^ m = 0, - ^m fc (-l) m ^_ m = 0. 

m 

When coupled with 

= V2 

and the orthonormality constraints, 

h[hi_2n = SnO 

I 

we get a system of equations that can be solved for the Daubechies-i^ scaling 
coefficients. The cases K = 1,2,3 have analytic solutions. These solutions 
are given in Table 1. 



Scaling coefficients for other values of K are tabulated in the literature [1]. 
With the exception of the Haar case (K = 1), there are two solutions which 
are related by reversing the order of the coefficients. 
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Given the scaling coefficients, hi, it is possible to use them to compute 
the the scaling function. While the Fourier transform method can be used to 
compute the Haar functions exactly, it is more difficult to use in the general 
case. 

An alternative is to compute the scaling function exactly on a dense set 
of dyadic points. This construction starts from the scaling equation in the 
form: 

(j,(x) = ^2V2hi(l>(2x-l). (105) 
i 

Let x = n to get relations between the values of the scaling function at 
integer points 

(j)(n) = J^v / 2^0(2n-/). (106) 
i 

Set m = 2n — I to get 

<t>{n) = ^h 2 „(t>{m) (107) 

m 

This gives the equations 

0(n) = ^#„ m 0(m) (108) 

m 

for the non-zero (j>(n) corresponding to n — 1, • • • , 2K — 2 where 

H nm = V2h 2n . m . (109) 

Eigenvectors of the matrix H mn with eigenvalue 1 are solutions of the scaling 
function at integer points - up to normalization. 

Rather than solve the eigenvalue problem, one of the equations can be 
replaced by the condition 

£>(n) = l (110) 

n 

which follows from the assumption that J ip(x)dx = 0. (The proof of this 
statement uses the fact that the translates of the scaling function on a fixed 
scale and the wavelets on all smaller scales is a basis for square integrable 
functions. Since 1 is locally orthogonal to all of the wavelets by assumption, 1 
can be expressed as a linear combination of translates of the scaling function. 
The normalization condition gives the coefficients of the expansion above.) 
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The support condition implies that only a finite number of the 0(n) are 
non-zero. This condition is independent of the orthonormality condition. 
For the case of the K = 2 Daubechies wavelets these equations are 

0(0) = >/2M(0) 

0(1) = v/2(M(2) + M(l) + M(0)) 

0(2) = ^(^10(3) + /i 2 0(2) + M(l)) 

0(3) = ^30(3) 

1 = 0(0) + 0(1) + 0(2) + 0(3). 

The first and fourth equation give 0(0) = 0(3) = (or h — hi — l/y/2 
which is the Haar solution). This also follows from the continuity of the 
wavelets, since and 3 are the boundaries of the support. The second and 
third equations are eigenvalue equations 

/ 0(1) \ _ / V2hi V2h \ (</>(!) \ 

{ 0(2) ; ^ V2h 2 V2h 3 ) { 0(2) ) ■ 

Instead of solving the eigenvalue problem for an eigenvector with eigenvalue 
1, use 

0(1) + 0(2) = 1 (112) 

with 

0(1) = >/2(M(2) + M(i)) 

to get 

0(l) = v / 2> o (l-0(l)) + M(l)) 
which can be solved for 

*>-TT^ho (113) 

and 

^ = l + y/Kho-h!) (lM) 

This gives exact values of the scaling function at integer points in terms of 
the scaling coefficients. This solution satisfies J2 n 4>(n) = 1. In this case 
there are only two non-zero terms. 
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In order to construct the scaling function at an arbitrary point x the first 
step is to make a dyadic approximation to x. Let m be an integer that defines 
a dyadic resolution. This means that we want the dyadic approximation to 
satisfy the inequality \x — x approx \ < 2~ m . For any m it is possible to find an 
integer n such that 

— <a;<-^. 115 

Writing this as 

n<2 m x <n+l (116) 

immediately gives 

n ■= [2 m x] = floor(2 m x) (117) 

where [] means greatest integer < 2 m x. 

Since the scaling function is continuous, for any e > we can find a large 
enough m so 

\<t>{x)- 0(^)| <e. 

In what follows we evaluate 0(n/2 m ) exactly Let x = n/2 m . We also assume 
that < n < 2K — 1 x 2 m , otherwise <p{x) = by the support condition (in 
this example we consider the case K = 2) . In order to evaluate <f>(x) note 
that the scaling equation gives: 

*(*> = * (£ ) = ^* = 

E ^^*(^r) = E ^ (£r - '0 

= E^'i*( n 2^/* ) ( 118 ) 

Repeating this process a second time gives 

(„ _ nm-l] _ om-2/ \ 
2 I 2 2 • (H9) 

Using the scaling equation m times gives 
(f>{x)= ^ m/2 h h h h ---h m Hn-2 m - 1 h-2 m ~ 2 l 2 ----2l m _ 1 -l m ). (120) 
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In this case the last expression is evaluated at integer values which gives (for 
the Daubechies K = 2 case): 

(f)(x) = 

^2 c h c h ---c m x (121) 



5 n _ 2 m-l/ 1 _ 2 m-2 i2 _...2i m _ 1 -i mi l 



V2hp 
1 + V2(h - h] 

1 - y/2tn 



+ 



1 + V2(h - h) 



(122) 
(123) 



where Ck ■= V%hk- 

This method generalizes to any value of K and any choice of scaling 
coefficients, hi. 

The scaling function and mother wavelet for the Daubechies wavelet are 
pictured in Figure 1. 
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Daubechies-3 Basis Functions 



Scaling Function 




1 2 3 4 5 

Figure 1. 



4 Daubechies Wavelets 

The Daubechies wavelets have two special properties. The first is that there 
are a finite number of non-zero scaling coefficients, hi, which means that 
the scaling functions and wavelets have compact support. The order-if 
Daubechies scaling equation has 2K non-zero scaling coefficients, and the 
support of the scaling function and mother wavelet function is on the inter- 
val [0, 2K — 1]. The second property of the order- K Daubechies wavelets is 
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that the first K — 1 moments of the wavelets are zero. 

The second property of the Daubechies wavelets is what makes them 
useful as basis functions. The expansion of a function f(x) in a wavelet basis 
has the form 

f( X ) = ^/mnt(l) fmn ■= / f (x)tfj mn (x)dx . 

mn 

If f{x) can be well-approximated by a low-degree polynomial on the support 
of i/j mn (x), then the vanishing of the low-order moments of ip mn {x) means 
that the expansion coefficient f mn will be small. On the other hand, as we 
will show in this section, the scaling function basis can be used to make 
local pointwise representation of low-degree polynomials. Since the scaling 
function basis on V m is equivalent to the wavelet basis on all scales, k > 
m, this means that the wavelet basis provides an efficient representation 
of functions that can be accurately approximated by local polynomials on 
different scales. For integral equations with smooth kernels, this means that 
the matrix representation of the kernels in a wavelet basis will be represented 
by a sparse matrix. 

The constraint on the moments of the Daubechies wavelets, 

'-0-^-1. (124) 
has important consequences. Eq. (j!24|) implies 



ip Qm (x)x l dx = J if)(x — m)x l dx = J ip(y)(y + m) l dy 

m l ~ k J ij(y)y k dy = I = • • • K - 1, (125) 



k\(l - k)\ 

which means that first K — 1 moments of the unit translates of the mother 
wavelet function vanish. Similarly, changing scale gives 

ipi${x)x l = J Dip{x)x l dx = —j= J ip{x/2)x l dx 



2 m/2 j 4){y)y l dy = l = 0---K-l. 



(126) 
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It straight forward to proceed inductively to show for all m and n that 

J^ nm (x)x l =0 l = 0---K-l. (127) 

This means that every Daubechies wavelet basis function is orthogonal to all 
polynomials of degree less than K, where K is the order of the wavelet basis. 
For the orthonormal basis of L 2 (M) consisting of 

{T n 0(x), D m T n tfj(x) : m < 0} (128) 

the only basis functions with non-zero moments with I < K are the scaling 
basis functions 

yW^O l = 0-K-l. (129) 

Although polynomials are not square integrable, we can multiply a polyno- 
mial by a box function b(x) which is 1 between X- and x + and zero elsewhere. 
The product of the box function and the polynomial is square integrable and 
is equal to the polynomial on the interval [ 1. It follows that 

p(x)b(x) = ^ C mn1pmn(x) = d n (j) kn (x) + ^ C mn lp mn (x) (130) 

mn n n m<k 

where p(x) is a polynomial of degree less than K and 

Cmn = / i) m n{x)p{x)dx (131) 



X- 



d n = <p n (x)p(x)dx. (132) 

J X- 

The moment condition means that the coefficients c mn = whenever the 
support of the wavelet is completely contained inside of the interval x + ]. 
Thus in the first expression the non-zero coefficients arise from end-point 
contributions and from many small contributions from wavelets with support 
that are much larger than the box. 

If k is set to correspond to a sufficiently fine scale, so the support of 
all of the wavelets is much smaller than the support of the box, then the 
second sum in (J13U)) has no wavelets with support larger than the width of 
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the box. The endpoint contributions only affect the answer within a distance 
A, equal to the width of the support of the scaling basis function, from the 
endpoints of the box. Inside this distance the only nonzero coefficient are 
due to the translates of the scaling functions. There are a finite number of 
these coefficients, and in this region they provide an exact representation of 
the polynomial. Specifically let 



I(x) = b(x)p(x) - 2^ d n (j) kn (x) + 2^ 2^ c mn ijj mn (x) (133) 

n n m<k 

then we have 

X— +A rx+ 

I(x) 2 dx + I I(x) 2 dx+ 

Jx + -A 
rx+ — A 

/ \ P (X) -J2 d n0kn(x)\ 2 dx. (134) 

Jx.+A 

Since all three terms are non-negative we conclude that 

rx+ — A 

/ \p( x ) ~ ^2d n (p kn (x)\ 2 dx = 0. (135) 

Jx.+A 

Since A is fixed by the choice of the support of the scaling function and x± 
is arbitrary we have 

/ \p(x) - ^2d n c/) kn (x)\ 2 dx = (136) 
Ja n 

for any interval [a, b}. Since p(x) and 4>(x) are continuous (we did not prove 
this for 4>(x)) and the sum of translates has a finite number of non-zero terms, 
it follows that 

P( x ) = ^2d kn (p n (x) (137) 

n 

pointwise on every finite interval. Since the box support is arbitrary this holds 
for any k. This establishes the desired result, that polynomials of degree less 
than K can be represented exactly by the finite linear combinations of the 
scaling functions (j) n {x). Since both bases in (jl30j) are equivalent, it follows 
that local polynomials can also be represented exactly in the wavelet basis. 

Figure 2. shows integer translates of the Daubechies 2 scaling function. 
Note how the sum of the non zero wavelets at any point in identically one, 
in spite of the complex fractal nature of each individual scaling function. 
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Daubechies-2 Scaling Functions Translated 




1 2 3 4 5 

Figure 2. 

Figure 2. shows a local representation of a constant function in terms of 
ling fucntions. 
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Figure 3. shows a local representation of a linear function in terms of 
scaling functions, while Figure 4. shows a local representation of a linear 
function. 
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Summing Daubechies-2 Wavelets to Represent a Linear Function 




-3 -2 -1 1 2 3 4 5 6 

Figure 4. 



Note that expansion in the wavelet basis gives all coefficients zero. This 
is not a contradiction because none of the polynomials are square integrable. 
The key point is that once one puts a box around a function, wavelets with 
very large support (large m) lead to many small contributions. 
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5 Moments and Quadrature Rules 



One of the most important properties of the scaling equation is that it can be 
used to generate linear relations between integrals involving different scaling- 
function or wavelet basis elements. In this section we show how the scaling 
equation can be used to obtain exact expressions for moments of the scaling 
function and wavelet basis elements as functions of the scaling coefficients. 
These can be used to develop quadrature rules that can be applied to linear 
integral equations. In section 6 the scaling equation is used to obtain exact 
expressions for the inner products of the these functions and their derivatives, 
which are important for applications to differential equations. We also show 
that these same methods can be used to compute integrals of scaling functions 
and wavelets over the different types of integrable singularities encountered 
in singular integral equations. 

Moments of the scaling function and mother wavelet function are defined 

by 

< x m >^= J <\>{x)x m dx < x m >^= J i){x)x m dx. (138) 

Normally these are integrated over the real line. For compactly supported 
wavelets this is equivalent to integrating over the support of the wavelet. 

A polynomial quadrature rule is a collection of N points {xi} and 
weights {wt} with the property 



r N 

< >., / o(. r ).<■"■ d.r = J2 x T w i ( 13!)) 



which hold for < m < 2N — 1. By linearity this means that 

N 

/ <p(x)P(x)dx = P{xi)wi (140) 
J i=i 

is exact for all polynomials of degree up to 2N — 1. 

In order to construct a quadrature rule we need to first compute the 
moments, and from these we can compute the points and weights. The 
moments can be constructed recursively from the normalization condition 

< x° >^= (x°, <P) = J dx(f)(x) = 1 (141) 
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using the scaling equation 



< x m >tf> = ( x m ,(f>) = {Dx m ,D<f)) 

11 m i 

1 1 ml 



r -k < x k 



Using ^ /i/ = a/2, and moving the k = m term to the left side of the above 
equation gives the recursion relation: 

m-l , /2JC-1 \ 

< x m > ,= L V" m ' y ^r- fe < x k >s . (142) 

Note that the right hand side of this equation involves moments with k < m. 
Similarly the moments of the mother wavelet function are expressed in terms 
of the moments of the scaling functions using eq. (|95p 



m . /2K-1 \ 

<x m >,= — y - m - - y M <x k >^. (143) 



2™^ k\(m-k)\ V2 

Since the scaling equation for the mother wavelet function relates the mother 
wavelet function to the scaling function there is no need to take the k = m 
term to the left of the equation; it is known from the first recursion. 

Equations (|141|) . (|142jl . and ()143j) give a recursive method for generating 
all non-negative moments of the scaling and mother wavelet function from 
the normalization integral of the scaling function. 

The moments for (pki = D k T l <p and ipki = D k T l ijj can be computed from 
the moments ()142j) and ()143|) using the unitarity of the D and T operators 



< x m > 4>kl = (x m ,D k T l 0) = {T~ l D~ k x T 



m . 

ml 



2 fc(m+l/2) ln - jm-n < % n > 

^— ' n\(m — n)\ 

n=0 v ' 
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and 



m . 

ml 



= 2 fc(m+l/2) V- ^ ,m-n < x « > 

^—^ n!(m — n) 

n=0 v 7 

Thus, all moments of translates and scale transforms of both the mother 
wavelet and scaling functions can be computed exactly in terms of the scaling 
coefficients. 

Partial Moments: Partial moments of the form 



pn 

< X m >0[ O :„]= / (p(x)x m dx 

Jo 



and 

-2K-1 



< X m > <t ,[ n ..2K-l]= / <j)(x)x m dx =< X m > (f) - < X m > < f >[ 0:n] 

J n 

for n G {1, • • • , 2K — 2} are also needed of numerical applications. 

These are can be calculated recursively in terms of the full moments using 
the scaling equation. First consider the order m = partial moments. Use 
the scaling equation 



^(x) = ^(|) = X)M(a:-0, 



which gives 

<j>(x) =^V2h l <j>(2x-l), 
i 

in the definition of the m = partial moment to obtain: 



< X° >0[Q:r 



pn pn 

/ (f)(x)dx = j <K 2;r ~ l ) dx - 

Jo j Jo 



Substituting y = 2x — I gives 

< x° > 0[O: „]= J2~^J t = J2^ <X ° >4>[-i,2n-i] ■ 
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The support condition of the scaling function (f>(x) implies that the lower 
limit of all of the integrals can be taken as zero which gives the relations: 



2n 



„0 



< >4>{0:n}— ^ C 2n -k < >4>[0:k] 
k=2n-2K+l 



where 



hj_ 
V2 



are non-zero for / = 0, • • • , 2K — 1. These equations are a linear system for 
the non-trivial partial 0-moments, m n =< x° > < f ) p. n }i m terms of the full 
0-moments < x° ><p= 1. These equations have the form 



M mn m n = v r , 



(144) 



where n, m : 1 — > 2fT — 2 and 



M mn S mn Cmn 



with 



and 



C mn = C 2m -n TR < K - 1, U = 1, • • • , 2m 

C mn = m < K — l,n = 2m + 1 , • • • , 2K — 2 

Cmn — C2m-n m — K — 1, K 

C mn = c 2m -n m> K,n = 2m-2K +1,---,2K -2 
C mn = m > K,n = 1, • • • ,2m - 2K 

v n = n < K 

2(n-K)+l 

Vn = ^ c k K <n <2K -2. 

k=0 

The matrix M := I — C has the form 

/ 1 - ci -c • • • 
-c 3 1 - c 2 -ci -c • ■ • 



I-C:= 



\ 



— C2K-3 —C 2 K-4 

— C2K-1 —C-2K-2 




-c 
-c 2 

"C2K-1 1 — C 2 i^_2 / 
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the vector m n , of partial moments, has the form 



m : = 



( <X° > ,£[0:1] \ 
< x° > 0[O:2 ] 



< X ><f>[0:K-l] 
< X° >4,[ :K] 



\ <X° > (f) [ :2K-2] J 

and the driving term has the form 



/ 



v : = 








c + Ci 



\ 



\ C H C 2 K-3 ) 

The solution of this system of linear equations gives the partial moments of 
order zero: 

m n =< X° >0[ O :n]= (I - C^v. 

Complementary partial moments are given by: 

< X° > < j > [n:2K+l]— l~ < X >0[O,n] • 

Higher order partial moments can be constructed similarly 



< or > 



n] := / (f)(x)x k dx 
Jo 

= / <K 2;r - l)x k dx 

;-n JO 



2K-1 



1=0 
2K-1 



£ 

1=0 



hi [ 2n ~ l 



<p(x)(x + lfdx 
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2A'-1 , k r 2n-l 

= V V - l k ~ m / <£fxk m cte 

(=0 m=0 

Rearranging indices, putting terms with the partial moments of the highest 
power on the left gives the following equation for the order k partial moments: 

2K~2 

^ (<W ~ -jk C rnr) < X k >0[ O: r]= W m (145) 
r=l 

where the inhomogeneous term w n in ()145|) is 



r=2n 

w n = 6(n>K) ^<x k >- 

2K-1 fc-1 



2fc 

2n-2K+l 



1=1 m=0 v ' 

which can be expressed in terms of the full moments of order k and partial 
moments of order less than k. The desired order-/c partial moments are 
obtained by inverting the matrix 



< x k > t 



Note that C matrix is identical to the C matrix that appears in equation for 
the 0-order partial moments. 

Having solved for the partial moments of the scaling function, it is possible 
to find the partial moments for the mother wavelet function < x m >^ k [o: n ] 
using 



1 "ill"' 

k\(m — k)\ 



< X m >^ k [o-nY— .,,„ . ; 9l X! HU _ Ml " X "' >( '.' (,: ~"-'- 



where the < x m ><f> k [o:2n-i] vanish for 2n — I < 0, are partial moments for 
< 2n — I < 2K — 1, and are full moments for 2n — l> 2K—1. This equation 
expresses the partial moments of the mother wavelet function directly terms 
of moments and partial moments of the scaling function. 

Given the moments and partial moments of 4>(x) and ip{x) we can solve 
for the partial moments of <p mn and tp mn in terms of the partial moments of 
(f>(x) and ip(x) by rescaling and translation. 
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Quadrature Rules: Given exact the expression for the moments it is pos- 
sible to formulate quadrature rules for integrating the product smooth func- 
tions times wavelet or scaling basis functions. The simplest quadrature rule is 
the one-point quadrature rule. To understand this rule consider integrals 
of the form 

J f{x)<t>{x)dx. (147) 
The quadrature point is defined as the first moment of the scaling function 



x\ := J x(f)(x)dx = fii. (148) 
With this definition we have 

/( a + ta)«*),l» = « + = « + *»,. (149) 
For orthogonal wavelets note that 

k m '■= / x(j)(x)4>(x — m)dx — / (y + m)(f)(y + m)(j)(y)dy 



x(f>(x + m)(f)(x)dx = k- m . 
It follows that 

^2mk m = ^2(-m)k m = 0. (150) 

m rrt 

For the Daubechies order K > 1 scaling functions 

m(j)(x — m) = x — ^ (151) 

m 

which gives 

O^J^*-^^/^-*^-^ (152) 

m ' m 

This means that fi2 = \A or 

J (j)(x)(a + bx + cx 2 )dx = a + bfii + cfi 2 = a + bx\ + cx\. (153) 
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This gives a one-point quadrature rule with point X\ = ji\ and weight W\ = 
1, that integrates the product of the scaling function and polynomials of 
degree 2 exactly. This is very useful and simple when used with scaling basis 
functions that have small support. 

More generally, given a collection of 2N moments of the scaling function 
or mother wavelet we can construct quadrature points and weights using the 
following method [?]. If {xi} are the (unknown) quadrature points define the 
polynomial 

N N 

P(x)=H(x-X l )=Y,PnX n 

i=l n=0 

where p^ = 1 and the other p n 's are unknown. Define the polynomials of 
degree n + m 

Q m (x) = x m P(x) 

for m = 1,---N — 1. By construction, for each m and Xi , Q m (xi) = 
because P(xi) = 0. 

If we require that the points {xi} and weights {wi} exactly reproduce 2N 
moments then it follows that 

N 

/ <p{x)Q m (x)dx = ^2 Qrn{Xi)Wi = (154) 

J 1=1 

because Q m (xi) = 0. The condition that the weights reproduce the moments 
give the conditions 



f N 

/ <p{x)Q m (x)dx = J2Pn< X n+m > , 

J 71=0 



and this must be equal to zero for each value m from m = to m = N — 1. 
This gives iV linear equations for the iV unknowns po ■ ■ -Pn-i' 



N 



Y,Pn < x n+m >4> = m = 1 • • -N; p N = 1 



n=0 



or 



N-l 



< * n+m ><pPn = ~< x N+m >^Pn m = 1 • • -N; p N = 1. 



n=0 
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Solving this linear system for the coefficients p n , using p N — 1, gives the 
polynomial P(x). 

Given the polynomial P(x) the next step is to find its zeros. The N zeros 
of the polynomial P(x) are the quadrature points Xi. Given the quadrature 
points, the weights are determined from the remaining N moments by solving 
the linear system 



for the weights, Wi. 

This shows how to construct the quadrature points and weights from the 
moments. In applications the linear equations for the coefficients p n are real . 
It follows that the zeros of P(x) are either real or come in complex conjugate 
pairs. In general it is desirable that the points are real and lie in the support 
of the scaling function. When this fails to occur it is best to simply assign 
real quadrature points that lie on the support of the scaling function. In 
doing this some accuracy is sacrificed, but it is easy to go to a higher order. 
Generally quadrature rules are used to integrate over the support of a scaling 
function of a small scale; normally a small number of quadrature points and 
weights is sufficient. 

For quadrature rules on a half-interval the partial moments, < x m ><^[o : oo]> 
need to be used near to generate a quadrature rule. 

Quadrature points are normally needed for different scaling basis func- 
tions, 4> mn {x). Points and weights for integrating <p mn (x) can be generated 
by scaling and translation. To see this consider a set of points and weights 
{xi,Wi} that satisfy 



N 





i=l 




It follows that 



(x m ,<j) nk ) = {x m ,D n T k <f)) = 2 n(m+1/2) (:r m ,T fc 0,) 
- 2 n{m+1 l 2 \(x + k) m ,<p) = Y,2 nm+n/2 Mxi + k)" 
= J2(2 n/2 m)(2 n (x l + k)) m . 



rn 
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If we define the transformed points and weights by 



w' l = 2 n ' 2 w l x' l = 2 n (x l + k), 

we get 

The new points and weights involve simple transformations of the original 
points and weights. 

While it is possible to formulate quadrature rules for both the wavelet 
basis and scaling function bases, it makes more sense to develop the quadra- 
ture rules for the scaling function on a sufficiently fine scale. This is because 
the scaling basis functions have small support, which means that the quadra- 
ture rule will be accurate for functions that can be accurately represented by 
low-degree polynomials on the support of the scaling function. 

Integral Equations: To use the quadrature rules to solve linear integral 
equations first consider the non-singular integral equation 



Let 



f(x) = g(x) + J K(x,y)f(y)dy. 



where (f) sn (x) are translates of the scaling function on a sufficiently fine scale 
s. Inserting the approximate solution in the integral equation gives 

n n 

Using the orthonormality of the <p S m{x) for different m values and a suitable 
quadrature rule gives the equation for the coefficients f m : 

fm = ^2 9{xim)wi m + ^ ^ Wl m K(xi m , X kn )w kn f n 
I n l : k 



or 



E 



~ S T J Wi m K(xi m , X kn )W kn 
l,k 



fn = ^2g(xi m )wi m . (155) 
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Note that no integrals need to be evaluated, except using the local quadra- 
ture rules. In addition the points and weights only have to be calculated 
for the scaling function on one scale - the rest can be obtained by simple 
transformations. 

While the scaling function basis on the approximation space V s is useful 
for deriving the matrix equation above, it is useful to use the multiresolution 
analysis to express the approximation space as 



The representation on the right has a natural basis consisting of scaling basis 
functions 4> s+r>m (x) on the larger scale, s + r > s, and wavelet basis functions 
ipmn{x) on intermediate scales s < m < s + r. These two bases are related by 
a real orthogonal transformation called the wavelet transform. Normally 
the wavelet transform is an infinite matrix. In applications it is replaced by 
finite orthogonal transformation that uses some conventions about how to 
treat the boundary. 

To solve the integral equation the last step is to use the wavelet transform 
on the indices m, n. This should give a sparse linear system that can be used 
to solve for f n . While the precise form of the sparse matrix will depend 
on how the boundary terms are treated, the solution in the space V m is 
independent of the choice of orthogonal transformation used to get the sparse- 
matrix form of the equations. 

Given the solution /„ it is possible to construct f(x) for any x using the 
interpolation 



This method has the feature that the solution can be obtained without ever 
evaluating a wavelet or scaling function. The wavelet nature of this calcu- 
lation appears in the quadrature points and weights, which are functions of 
the scaling coefficients. 

Scattering integral equations have two complications. First the integral 
is over a half line. Second, the kernel has a l/(x ± it) singularity. 

The endpoint near x = of the half line can be treated using special 
quadratures for the functions on the half interval. If there is a <p n with 
support containing an endpoint, the 5 mn in (|155|) needs to be replaced by 



V s = Ws+l © Ws+2 © ■ ■ • © W s+r © V s+r . 




n 



k 




43 



which is not a Kronecker delta when the support of m and n contain 
the origin. These integrals can be evaluated using the same methods that 
were used to calculate moments on the half interval. We use the scaling 
equations and the ort honor mality when the support of both terms are in the 
half interval. Specifically for a and b integers 



N$= f O^oMHr 

J a 

<(x)(f>(x + i — j)dx 

/O—l 
(f)(2x - l)(f)(2x + 2i- 2j - l')2dx 

i-2(b-i)-l 

^2 h ih' \ 4>(x)(j)(x + 2i-2j + l- l')dx 

i J 2(a—i) — l 



f 

J a- 



E, , M 2a-2i-l,2b-2i-l) 
nini>M ,-2i+2j-i+i' ■ 



i,v 

When either function has support inside the interval this is a Kronecker delta. 
These equations are linear equations that relate these known elements to the 
unknown elements where the support overlaps an upper or lower endpoint. 
These formulas simplify if one of the endpoints satisfy a = oo or b = —oo. 
The final relations are 

iKja—i,b—i \ v 7 ii i\j-2a—2i—l,2b—2i—l 

iV 0,i-i - tl l tl l I ^0-2i+2j~l+l> ■ 

1,1' 

Note that N , k = 1 if k = 0, N 0)k = if k > or k < -{2K - 1). It is 
non-trivial for —(2K — 2)<k<0. This gives a linear system for the overlap 
coefficients, N i:j . 

For scaling functions that overlap x = the equation becomes: 



E 



No,n-m — Wi m K(xi m , X kn )W kr , 



l,k 



fn = ^2g(xi m )wi r 



where the xi m ,wi m indicate that for m satisfying 2K — 2 < m < the 
quadrature points and weights need to be replaced by the ones for the half 
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interval. In this case overlap matrix elements N mn need to be computed on 
the scale dictated by the approximation space 14 . 

Mapping techniques are valuable for transforming the equation to an 
equivalent equation on a finite interval and for treating singularities. For 
example, the mapping 

by - a 

X = -X -- 

a o — y 

transforms the domain of integration to [a, b] and a singularity at x = x$ to 
the origin. What remains is a mechanism for treating an integrable singu- 
larity. The first step is to use a mapping, like the one above, to place the 
singularity at the origin. After mapping, the relevant integrals for a l/(x—x ) 
singularity, when using the subtraction technique discussed below, are 

Un):= j£^m dx . 

Using unitarity of D gives 

I m (n) := [ D(D m T n <P(x))D-dx = 
J x 

2 r D m+1 T n (j)(x) 



I 



-dx 



y/2 J x 

2K-1 

>/2/i,/ ro (2n + Z). 



1=0 



The equations 



2K-1 

Im{n) = ^hlm&n + l) 

1=0 

give linear equations relating the integrals with singularities to the integrals 
with no singularities. The singular terms for the order-K Daubechies scaling 
functions are 

mO (-l), mO (-2), • • • (p m o(-2K + 2) 

The endpoint terms, m o(O) and <p m o(—2K + 1) are not singular because 
4> m (x) must be continuous at the endpoints. 
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We found that these equations are ill-conditioned, but they can be sup- 
plemented by 

o = ^v J k ^ mn ^ 

which has the form 



= J2 V ^^-dx, (156) 
J-k x 

= (157) 



where the integrals I^(n) include partial integrals when n is such that the 
support of mn (x) contains the points k or —k. These linear relations relate 
the singular integrals to the non-singular ones. For 4> mn {x) with support far 
enough away from the origin and not containing k or — k the integrals can 
be expressed in terms of the moments: 

^T^-dx = / T n (f)(x)D- m -dx = 2~ / (j){x)T- n -dx 



1 —dx = 2-f 1 -jr(- 1 -) l (x\ 



For large values of n this series converges rapidly Similar methods can be 
used for values of n where k or —k is in the support of (f> mn (x). In this case 
the full moments need to be replaced by the appropriate partial moments. 
For singularities of the form 1/ (x ± ie) equation ()156|) is replaced by 

dx 



m 



x±iO + 



Using this with the wavelet expansion 

n 

provides the needed additional equation, 

n 

which replaces (|157|) . The result of solving these linear equations is accurate 
approximations to the integrals 

(158) 
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To use these to solve the integral equation consider the case m = 0: 



<Pn{y)dy 

y 

^ K(x, xi) - K(x, 0) n\r ( \ 

= w i + K K X , 0) J o(«) 

i Xl 

In applications the Io(m) should be computed for m values far from the 
singularity using the series expansion. The equations relating the Io(m) are 
used to calculate the remaining 7 (^)' s - 

Similar methods can be used to treat other types of integrable singulari- 
ties. For example, for logarithmic singularities define 

I (n) := J (f) n (x)\n(x)dx 

The scaling equation gives the linear relations 

I (n) = (T n (PM) = (DT n (P,D\n) 

= _^ [(T 2n jD0!ln) _ ln(2)(T 2n jD0)1)] 

= 7f ^>/o(2n + 0-ln(2) 

In this case, because the singularity is integrable and the value of the inte- 
gral is unambiguous, we do not need an additional equation to specify the 
treatment of the singularity; however the function is multiply valued, so the 
computation of the input integrals far from the singularity should reflect the 
choice of Riemann sheet. 

The linear equations above relate the integrals Io(n) that overlap the 
singularity to integrals far away from the singularity, which may or may not 
have support containing the endpoints ±a. These terms serve as input to the 
linear system and can be computed with the moments and partial moments 
using expansion methods. For the case that the support of <f> n (x) does not 
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Table 2: Singular Integrals 



K = 


2 






n = 


-2 


f(f)(x 


— n) In 


x\dx 


0.456927033732831 


n = 


-1 


f(f>(x 


— n) In 


x\dx 


-1.64215549088219 


K = 


3 






n = 


-4 




— n) In 


x\dx 


1.15737952417967 


n = 


-3 


f<f>(x 


— n) In 




0.750468355278047 


n = 


-2 


f<f>(x 


— n) In 


x\dx 


0.315624303943019 


n = 


-1 


f(f)(x 


— n) In 


x\dx 


-1.83646456399118 



contain ±a the following expansion can be used when the support of (f) n (x) 
contains positive values of x: 

h{n) = J <pn{x) \n(x)dx = J <p(y) ln(n(l + y/n))dy 

= ln(n) - Y ( ~ ir < ^ >4> . (159) 

m=l 

This expansion converges rapidly for large n. When |n| is large with n < 
we use 

ln(n) = ln\n\ + i(2m + l)n 

and m is fixed by the treatment of the integral near the origin. In the 
case that one of the endpoints ±a are contained in the support of <p n the 
moments, a similar expression can be used to approximate the integrals Iq{ti), 
except the moments, including the 0-th moment multiplying ln(n), need to 
be replaced by partial moments. 

For the case of the Daubechies K = 2 and K = 3 wavelets the solution 
of these equations in given in Table 2. The numbers in the table are for 
f 4>n{x) In \x\dx for the 4> n (x) with support containing x — 0. 



6 Derivatives and Differential Equations 

The scaling equation can also be used as a tool to solve differential equations. 
Consider the following approximation of the function f(x) given by 

f(x)~Y,f n< f>™(x) ( 16 °) 

n 
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where (/) mn (x) = D m T n (p(x) is the scaling function basis on the approximation 
space V m . As m — > — oo this representation becomes exact. 

For the purpose solving differential equation we want to calculate approx- 
imate derivatives of f(x) on the same approximation space 

f'( x ) ^d' n (j) mn (x). 

n 

f"( x ) ^d„(f) mn (x). 

n 

The orthonormality of the scaling basis functions can be used to find the 
expansion coefficients d' n and d" n : 

d 'n '■= / 4>mn{x)f\x)dx = - <f)' mn (x) f{x)dx. (161) 



< := J <P m n(x)f"(x)dx = J 4>l n {x)f{x)dx. (162) 
Using the expansion (jl6()j) in p61j) gives 

d'n = - '^2(4>'mni 4>rnn>) fn> 



and 



The coefficients ((p' mn , 4> mn >) and (0^ n , <f> mn ') are needed to compute the ex- 
pansion coefficients d' n and d" n for the derivatives in terms of the expansion 
coefficients f n > of the function. 

Given these linear relations between the coefficients d' n and d", the so- 
lution of a linear differential equation can be reduced to solving a linear 
algebraic system for the coefficients /„. The size of this system can be re- 
duced by employing the wavelet transformation. The method of solution 
depends on the type of problem. Standard methods can be used to enforce 
boundary conditions; the only trick is that all of the basis functions with 
support that overlaps the boundary should be retained. 

The goal of this section is to show that the scaling equation can be used 
to compute the needed overlap integrals. 
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To proceed we first consider the simplest case where the scale m — 0. 
Using unitarity of the translation operator gives the following relations 

{4>' n ,M = {<t>'An'-n) = (<f>' n - n ' , 4>) 
(<f>l<f>n>) = (<f>",<f>n>-n) = (<f>'L n ',<f>). 

In addition these coefficients have the following symmetry relations 
= - y <t>\y)<t>(y + n)dy = -((/>', 4>- n ) 

and 

5 ((!/'. (!)- ri ) 



which follow by integration by parts. 

The overlap coefficients can be computed using the scaling equation and 
the derivatives of the scaling equation: 



2K-1 



4>(x) = V2 h t (j>{2x - I) 

1=0 
2K-1 

4>'{x) = 2V2 h4>\2x - i) 



1=0 
2K-1 

f(x) = (2 2 )V2^M>-0. 
1=0 

We first consider the computation of the coefficients (<f>' (x) , <f> n ) . 
For a n defined by 

a n := ((f)'(x),(j) n ) 

this leads to the following linear relations among the overlap coefficients 

2K-1 „ 

a n = 4 h i h v I ^ x - 2n - l)<j>'{2x - l')dx 

l,l'=0 J 

2K-1 „ 2K-1 

= 2 hhv <f>(y-2n-l + l')4>'(y)dy = 2 ^ h l h l ,a 2n+ i-i>- (163) 
1,1' =0 J 1,1' =0 
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Since both 4>(x) and 4>'(x) have support for < x < (2K — 1), the non-zero 
terms in the sum are constrained by 

-{2K - 1) < 2n + I - I' < 2K - 1. 

For the second derivative these equations are replaced by 

2K-1 

a n : = (0"(x),0 n ) = 8 ^ hih v a 2n +i-v (164) 

l,l'=0 

These linear equation are homogeneous and must be supplemented by a 
normalization condition. For the Daubechies wavelets of order K > 1 we 
have the expansion 

n n 

where the expansion coefficients are 

K. = / m*)** = » + / «* - »><* - »>* = »+ < * >* 

and 

\j' n = J 4> n (x)x 2 dx = J <f)(x) (x + n) 2 dx = n 2 + 2n<x> ( f > + <x 2 > ( f, 
= n 2 + 2n < x > + < x >\= (n+ < x > lj> ) 2 . 

Thus 

X =< X >0 + ^2 n 4>n{x) X 2 =< X >\ +2(x- < X ></,) < X > ( f l + ^2 n2< Pn{x) 

n n 

These equations can be differentiated to get 

1 = J>^(x) 2 = ^n 2 C(x) 

n n 

Multiplying by (f>(x) and integrating gives the desired inhomogeneous equa- 
tion 

n n 
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- ^Jn(0 n , <p') = - ^2na n (165) 

n n 

and 

2 = ]>> 2 < (166) 

n 

Equations ()163J) and ()165|) or (|163|) and ()166|) are linear systems that can be 
used to solve for the coefficients a' n and a!' n . 

In general, it is desirable to expand a function using a scaling basis associ- 
ated with a sufficiently small scale m, In addition, for efficiency it also useful 
to use the basis on the approximation space V m consisting or the scaling 
functions on the scale m + k and wavelet basis functions on scales between 
m + 1 and m + k. Finally, one needs to be able to treat higher derivatives. 
Generalizations of the methods can be used to find exact expression for all of 
these quantities expressed a solutions of linear equations involving the scaling 
coefficients. For the higher derivatives it is necessary to use a Daubechies 
wavelet of sufficiently high order. The number of derivative of the wavelet 
and scaling function basis increases with order. 

A necessary condition for the solution of the scaling equation to have k 
derivatives can be obtained by differentiating the scaling equation k times, 
which gives 

<P {k \x) = V22 k J2 h i<P k ( 2x - 1 ) 
i 

Letting x = m and n = 2m — I gives the equation 

^(m) = v / 22 fc ^/ i2m _ n /(n) 
i 

where the non-zero values of 4> k (n) satisfy < n < 2K — 1. For this equation 
of have a solution the matrix H mn := /i2m-n must have eigenvalue 2~ ( - k+1 / 2 \ 
This is a necessary condition for the basis to have k derivatives. When the 
k-th derivative exists, it can be computed up to normalization by iterating 
the Fourier transform of equation (jUJ). The method used to compute wavelets 
at dyadic points can also be used with the above equation to compute the 
A;-th derivatives of scaling functions and wavelets at dyadic points. 

The scaling equation can be used to exactly compute all of the expansion 
coefficients. In order to exhibit the key relations it is useful to use operators: 

Df(x) = (167) 
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Tf(x) = f(x - 1) (168) 

Af{x) = ^ {X) - (169) 
Direct computation shows the following intertwining relations 

AD = ^DA (170) 

DT = T 2 D (171) 

AT = TA (172) 

A f = -A T f = T- 1 D ] = D~\ (173) 
We also have the scaling equations: 

D ( j ) = J2hiT l <j ) (174) 
i 

Dij = J2giT l (p- (175) 
i 

Using the operator relations above give 

DA r 4> = 2 r ^2 h{T l A r (j) (176) 
i 

DA r ip = 2 r 9i Tl ^ r< P- (177) 
i 

The different expansion coefficients can be expressed in terms of these 
operators as 

(</w, A r <p mn ) = (D m 'T">, A r D m T n (fi) (178) 

(Vw, A r cf> mn ) = (D m 'T n '^, A r D m T n (f)) (179) 

(fa*, A> m „) = {D k 'T n '<t>, A r D m T n ip) (180) 

(Vw, A^ mn ) = (D k 'T n '4>, A r D m T n 4>). (181) 

The following steps are used to evaluate these coefficients: 

1. Move all of the factors of D to a single side of the equation. Choose 
the side where the power of D is positive. 

2. Move the D's through all derivatives. 

3. Use the scaling equations to eliminate all of the D's. 
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4. Move all of the T's to the left side of the scalar product. 

Using these steps all of the overlap coefficients can be expressed in terms 

of 

(0n,A» (182) 
ty>n,A» (183) 

(0n,A>,) (184) 

(V>n,A^) (185) 

The scaling equation can be used to express all of the ip terms in terms of 
the 4> terms. The result is at all of the coefficients can be expressed in terms 
of the coefficients 

(0n,A» (186) 

We have shown how to compute these for r = 1 and 2. The coefficients for 
larger values of r can be obtained by solving the system: 



2K-1 



a« := (0"(x), n ) = 2 2 - 1 hh v atl + i-v (187) 

l,l'=0 

7 Galerkin for Scattering 

We want to find the solution to the s-wave Schrodinger equation 
h 2 1 d 2 

—[rij(r)} + V(r)ij(r) = Eij(r) (188) 

2m r ar z 

for a particle with mass m and energy E = h 2 k 2 /2m, where ip(r) has the 
asymptotic form 

rip(r) — > sin(£r + 5) , (189) 

r— »oo 

and rif)(r) is zero at the origin. Equation ()188j) can be rewritten in the form 

1 d 2 

r^[r^(r)] + U(r)ip(r) = k 2 ^(r) , (190) 

r dr 2 
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where 

2/77, 

U(r) = —V(r) . (191) 

To solve Equation (jl90j) we choose a complete set of basis functions <p n (r) 
and write 

N 

<Kr) = J2 a n<Pn{r)- (192) 

71=1 

To solve for the expansion coefficients a n we first multiply Equation ()19()|) by 
'" 2 0m(^) and integrate from to R. This gives 

" R d 2 f R f R 

f 4>m(j)—;[r ip(r)]dr+ / (f) m (r) U(r) ijj(r)r 2 dr = k 2 / (f) m (r) ijj(r) r 2 dr . 



o dr 2 j j 

(193) 

Using integration by parts, the first term in Equation ()193j) can be written 

as 



R A2 







R r R d , ,j 



d 2 d 
-/ r(p m (r) — [r^(r)]dr = -r(t) m (r) — [r%l){r)}^ + j ^MmW]^M(r)](lr. 

(194) 

Now we set 

[rV>(r)] =1, (195) 



r=R 



which corresponds to a change in the normalization of ^(r), and use rip{r) 
is zero at r = to write 

f R d 2 f R d d 

(196) 

Thus, Equation ()193|) can be written in the form 

f R d d f R 

I —[r4> m (r)] — [rip(r)}dr+ / <p m (r) U(r) ip(r)r 2 dr 
o dr dr J 

r R 

- k 2 / m (r) ${r) r 2 dr = R <f> m (R) . (197) 
Jo 

Substituting the expansion for %jj{r) in Equation ()192|) into Equation ()197|) 
gives 

N f f R R 

S a? M / M r ^rn(r)]£[r (j) n (r)\ dr + f* <p m (r)U(r) <p n (r)r 2 dr 



n=l 







k 2 j*<P m (r)<p n (r)r 2 dr\ = R<f> m (R) , (198) 
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which can be written as the matrix equation 

N 



5> 

71=1 



ran 



Urn. • 



(199) 



Given the solution to Equation (j!99|) we need to determine the normaliza- 
tion and the phase shift. For the new boundary condition given in Equation 

ri/j(r) — > A sin(kr + 5) , 



Thus, we need an additional equation to use with Equation 
the integral 



R 



sin(£r) U(r) ip(r) rdr 

f R (Id 2 ) 

J S in(£r) )- — [rij(r)} + k 2 ij(r)jrdr 



d R 
sm(kr) — [r ip(r)] — kcos(kr)[r ip(r)] 

Rr/i d2 ^ 2^ . „ > 

-— - + k smlkr) 
r dr 2 / 



rip(r) dr 



= kA [sm(kR) cos(kR + 5) - cos(kR) sin(kR + 5)] 
= —kA sin S . 

From Equations (|195j) and (|2(J(J|) we get 

kA cos(kR + 5) = kA cos S cos(kR) — kA sin S sin(kR) = 1 
which using Equation (|2(J1|) can be written as 

c l-Ism(kR) 
kA cos o = t— - — . 

cos(A;ii) 

Finally, from Equations ()201|) and ()203|) we find 

. — I cos(kR) 

tan d = — - . 

1 - Isin(kR) 

Given 5 the value of A can be found using Equation (|2()lj) . 



(200) 
We use 



(201) 



(202) 



(203) 



(204) 
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Eckart Wave Function 
To test the method we use the Eckart potential 

V(r) ~ ^ WX2e ~ XT (205) 
V{r) ~ 2m + 1)2 ' ^5) 

which has the analytic solution 

[{Ak 2 + A 2 ) + {Ak 2 - \ 2 )/3e' Xr ] sin{kr + 5)- Ak\(3e~ Xr cos(Ax + 5) 
r ^ (r) = {Ak 2 + A 2 ) {(3e~ Xr + 1) ' 

where the phase shift is given by 



If (3 is chosen to have the value 



the potential will have a bound state with the energy 



The boundstate wave function is given by 



8 Wavelet Filters 



(206) 



5 = aictm| A) +arctan (|^lI 1 . (207) 



B B = -*£. (209) 



, . , 2V2^sinh(Ar/2)e- Kr . . 



The wavelet transform is the orthogonal mapping between the scaling func- 
tion basis on a fine scale and the equivalent basis consisting of scaling func- 
tions on a coarser scale and wavelets at all intermediate scales. The wavelet 
transform can be implemented by treating the scaling equation and the equa- 
tion defining the wavelets as linear combinations of scaling functions for a 
finer scale, as low and high pass filters. This has the advantage when most 
of the high-frequency information is unimportant. 
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The wavelet filter is defined using the scaling relations 

1 2fc-l 

D0(x) = - 7 =0(|) = g^rV(x) (2ii) 

and 

where 57 = (— l)'/i 2 k-i-j- Using Equations (|211|) and (|212|) we can write 

fx) = D J T m <j)(x) 

2fc-l 

= hiD j T m D- l T l (t){x) 



'3,m\ 



1=0 
2k-l 



hD j - l T 2m T l (j){x) (213) 

2fc-l 

^ hl(f}j- lt 2m+l{x) 



1=0 



and 



2k-l 



giD>T m D- l T l <P{x) 

1=0 
2fc-l 

^ gi D j - l T 2m T l (j){x) (214) 



Z=0 
2fc-l 



i=0 



where we have used 



TD- l ct){x) =TV2(p(2x) 

= V2(p(2x - 2) (215) 
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£/>,/, 

I J 



and 



to write T m D~ x = D~ l T 2m . 

We use the orthonormality of the functions <pj tm (x) and ipj tm (x) to obtain 
the inverse relation. Since Vj_i = Vj © Wj, we can write 

<t>j-l,m{x) = ^ a n<fij,n( X ) + ^j,n{ X ) > ( 216 ) 

?i n 

Using the scaling relations, we find 

a n = J (p j<n (x) (pj-l-^x) d,X 

<Pj-l,2n+l{ x ) (f)j^ m (x) dx 

= £><Wm ( 21? ) 

I 

= h m _2n (218) 

K = J i>j,n{x) (pj-l,m(x) dx 

= ^.91 (f)j-l,2n+l(x)(f) j - 1:m (x)dx 
I J 

2n+l,m 

(219) 

i 

= 9m-2n ■ (220) 

Thus, we find 

4>j-l,m{x) = ^ hm-2n(f)j,n( X ) + ^2 9rn-2n^j,n{ X ) ■ ( 221 ) 
n n 

An alternate derivation is to use the orthonormality relations 

hihi +2m = <5m0 (222) 

X] Ml+2m = 5m0 (223) 

X>0+2m = O, (224) 
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derived in the Wavelet Notes. Now using the scaling relations in Equation 
dUSD gives 

4>j-l, m {x) = a n hl(j)j-l t 2n+l{x) + &n ^ gi(j)j-i )2 n+l(x) ■ (225) 
n / n 2 

Taking the inner product with <f>j-x,k gives 



n I n I 

= 22 a nhk-2n + 22 bn 9k-2n ■ (226) 
n n 

Multiplying Equation (|226|) by /i m _2n' and summing over m gives 

hk-2ri — a ' n hk-2rihk-2n + &n hk-2n'9k-2n 

n k n k 

= Gl„<W (227) 



n' j 



where we have used the relations (|222|) and (|224|) . Multiplying Equation 
()226|) by gk-2n' and summing over gives 

9k-2ri = On 9k-2n'hk-2n + &n 9k-2n'9k-2n 

n k n k 

= J2 bn6 ™' ( 228 ) 

n 

= K' , 

where we have used the relations ()223|) and ()224j) . 

Given the expansion of f(x) in terms of the scaling functions 

2 J -1 

f{x) = Y,c,, n <p hn (x), (229) 

n=0 

we can use Equation (|221j) to write the expansion in the form 

2 J -1 2 J -1 

f( X ) = c j,n y hri-2m 0j+l,m(^) + ^""^ C j,n ^ ^ 9n-2m'4 ) j+l,m{ x ) 

n=0 m n=0 m 

2 J-i_ 1 2- 7 - 1 -! 

= 2j C j+ltm (j) j+ltm (x) + 22 dj+l,m^j+l,m(x) , (230) 
m=0 m=0 
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where 



and 



2m+2p-l 

9+1,™ ^ h n —2m Cj,n 

n=2m 

2m+2p-l 

^j+l,m ^ ^ 9n—2m Cj,n j 

n=2m 



(231) 



(232) 



where we use the periodic wrap-around condition c^j+j = c^. Equations 
(|231|) and (|232j) can be written as a matrix equation, which for J = 3 has 
the form 



/ 


C j+1,0 






( h 


hi 


h 2 


h 











\ 




I 


9,o 


\ 




9+1,1 












h 




/i 2 














9,i 






9+1,2 


















h 


hi 


h 2 


h 






9,2 






9+1,3 






h 2 


ka 














ho 


h 






9,3 






4/4-1,0 






9o 


9i 


92 




















9,4 






^i+i,i 












9o 


9i 


#2 


93 












9,5 






dj+1,2 




















9i 


92 


#3 






9,6 




\ 


dj+1,3 


) 




\ 93 


94 














9o 


9i J 




\ 


9,7 


/ 



(233) 



for the Daubechies p = 2 wavelets. Repeated application of the filter trans- 
form to the remaining Cj+i jm gives 



/ 


9,0 


\ 




/ 


9+1,0 






( 


9+2,0 


\ 




( 


9+3,0 


\ 




9,i 








9+i,i 








9+2,1 








4+3,0 






9,2 








9+1,2 








^i+2,0 








dj+2,0 






9,3 








9+1,3 








dj+2,1 








dj+2,1 






9,4 








dj+i,o 








"i+i,o 








dj+i,o 






9,5 








dj+i,i 








dj+i,i 








dj+i,i 






9,6 








dj+1,2 








dj+1,2 








dj+1,2 




V 


9,7 


/ 




V 


dj+1,3 


) 




\ 


dj+i,3 


J 




V 


dj+i,3 


I 



(234) 



The reverse transform can be obtained by substituting Equations ()213j) 
and ()214|) into Equation (j230J) . This gives 



2 J ~ 1 -l 



2fc-l 2 J ~ 1 -l 2k-l 

'j+l,m 

1=0 m=0 1=0 



2k-l 



m=0 



) j,2m+/(3 ; ) 
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h n -2m Cj+l,m 0J,n(^) + 9n-2m dj+l t m 4>j,n{x) 

n m n m 

^ ] 9> < / ) j.«(- r ) > 



(235) 



where 



9,™ ^ ^ fori— 2m Cj+l,m ^ ^ 9n—2m dj+l.rn 
m m 

For the example used in Equation ()233|) . this gives 



(236) 
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9,3 
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93 


9i 
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93 


9i 
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h 3 


hi 








93 


9i J 
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dj+i,3 
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which is the transpose of the matrix in Equation 



9 Wavelet Transfrom 

The wavelet transform is by its nature local at each level and therefore admits 
an implementation in which the data to be transformed can be placed in a 
buffer instead of storing the entire data set at once. This significantly reduces 
the amount of storage space required for applications involving compression. 

In the one-dimensional case, the J-level wavelet transform can be com- 
puted by buffering O(J) nonessential elements or the full transform can be 
computed buffering 0(log(iV)) elements. The standard form for the two- 
dimensional transform of an N x iV matrix can be performed by buffering 
only 0(N\og(N)) elements. In general, a D-dimensional wavelet transform 
can be computed by only storing 0(N D ~ l log(iV)) elements. This buffered 
wavelet transform can be applied to any type of data that can be input or 
computed in series. Some notable examples include the compression of time- 
series data and applications to solutions of integral equations. Below, we 
will explain the exact implementation of the transform including the buffer 
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and the extension of the method to two dimensions. The extension to ar- 
bitrary dimension is straightforward from the two dimensional case. First, 
we will layout the terminology that will be used throughout, we adopt the 
filter viewpoint since it makes the explanation of the buffering procedure 
more clear. But this is equivalent to the linear algebra viewpoint and we will 
attempt to explain the procedure in this language as well. 

From the filter viewpoint, the wavelet transform is a convolution of the 
data set and two vectors h and g followed by a decimation. This is equiva- 
lent to a convolution that proceeds by steps of two instead of one. For the 
Daubechies family of wavelets, both of these filters have a length L = 2K—1. 
The convolution with h produces what is called the father set and the con- 
volution with g produces the mother set. We denote the data set to be 
transformed as A and use brackets to denote subscripts. In all contexts be- 
low, the one-dimensional length of the data set will be iV = 2 n where n is an 
integer. That is the data runs from A[0] to A[N — 1]. The typical procedure 
to deal with a finite data set is to periodize the data over the boundary such 
that A[N] =A[0],A[N + 1] = A[l], ■ ■ ■ A[N + L - 3] = A[L - 3]. 

Now, if we write out the first level of the transform as a matrix we can 
see that it is banded with a bandwidth L corresponding to the convolution 
operation. The second level of the transform is identical to the first except 
that it only acts on the father set of data, i.e. the transform on the moth- 
ers is the identity. This corresponds to the fact that all of the information 
about coarser levels is contained in the father functions. The mother func- 
tions form an orthogonal subspace to the fathers and mothers on all higher 
levels. Using this knowledge we can immediately perform a thresholding pro- 
cedure on the mother set without affecting the rest of the data in any way. 
The father set can simultaneously be transformed and the resulting mothers 
thresholded as well. Iterating this procedure on all the relevant levels forms 
the basis for the buffered wavelet transform. Of course, one must have an 
a priori thresholding scheme to accomplish this. The simplest such example 
is an absolute threshold. In this scheme, one chooses an epsilon and all el- 
ements with a magnitude less than this epsilon. Other more sophisticated 
thresholding procedures exist as well, such as procedures based on the level 
on the transform. The important fact is that one cannot have a procedure 
that depends in any way on the final transformed data set. Examples of such 
procedures would be based on the relative size of the transformed elements or 
a threshold that keeps a certain number/percentage of the final coefficients. 
In many applications, the absolute thresholding is an acceptable method. 
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Now, we will explain the detailed implementation of the one-dimensional 
transform. One begins by computing the elements A[0] to A[L — 3] of the 
data set. As noted above, these elements are necessary for the periodic 
boundary conditions and form a boundary buffer that must be saved until 
the end of the calculation. Now, elements can be added to a moving buffer 
of length L that constitutes the heart of the procedure. After the elements 
A[L — 2] and A\L — 1] are computed and placed in the moving buffer, one 
can begin transforming the data set. Convolving this data set, including the 
boundary terms, with h produces the first member of the father set F[0] and 
convolving with g produces the first member of the mother set M[0\. As 
described previously, this mother element can be immediately thresholded 
and placed in the final output vector. The father element is considered the 
beginning of a new data set to be transformed and is placed in the boundary 
buffer corresponding to the next level of the transform. One then proceeds 
to compute two more elements and convolve A[2] to A[L + 1] with h and g. 
This produces F[l] and M[l], which are treated the same way as before. We 
continue in this manner, calculating more elements and convolving, until we 
have computed the element A[2L — 3]. The moving buffer is now full and we 
have reached the interior of the data set. When we compute the next element 
of the data set we can discard the last element of the moving buffer and shift 
all the elements of the buffer one place. The new element is then appended to 
the moving buffer. Discarding the last element is justified by the fact that all 
the information in that element is represented by the corresponding father 
and mother data sets due to the equivalence of the subspaces. The name 
moving buffer is clear since this buffer can be viewed as scanning the interior 
of the data set by moving over it. This process continues, the shifting and 
convolving, until the end of the data set is reached. When the end of the data 
set is reached we simply make the data set periodic using the boundary buffer. 
This process is simultaneously carried out at each level. Now counting the 
elements in each buffer we see that in each level we must store L — 2 elements 
in the boundary buffer and L elements in the moving buffer. So, for J levels 
we must store J{2L — 2) elements. This gives us our size of O(J) where the 
coefficient depends on the length L of the wavelet filter as is common with 
most wavelet algorithms. 

In many wavelet applications, the data vector to be transformed will be of 
length N = 2 n and a wavelet transform of level J = n will be computed. In 
this case, the number of elements stored in the buffers will be of 0(log(iV)). 
A minor point to note is that for wavelet filters of length > 2 the last few 
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levels will not be filled completely. As a programming point one can either fill 
the buffers periodically or just periodize the convolution. Both procedures 
are equivalent and consistent with the periodic wavelet transform. Also note 
that the number of operations has been increased by the shift operation, but 
is still of O(N) which is the case for the standard wavelet transform. The 
standard procedure to perform the wavelet transform on a two dimensional 
data set is to first transform the rows of the matrix and then transform the 
columns. Alternatively, one could transform the columns and then the rows. 
Both are equivalent as can be seen by writing out the transform as a matrix 
multiply and noting the associativity of matrix multiplication. 

To perform the buffered wavelet transform on a two-dimensional data 
set we calculate the data column by column. Each row has a separate set 
of buffers associated with it. We can view this as a strip that scans the 
matrix in much the same way as the moving buffer did in the one-dimensional 
case. Each of these buffers behaves in exactly the same way as the one- 
dimensional case, except the output is handled differently. The first output 
from the buffer associated with row one is placed in two vertical buffers. 
FS[0] and MB[0] the B stands for blank since these are internal buffers 
that have no outside significance. Both of these outputs must be saved 
because they contain information about the columns of the matrix. Row 
two then produces F.B[1] and MB[1], and so on continuing down the rows. 
The transform procedure is applied to the vertical buffers, which produce 
output FF, FM, MF, and MM. The output MM can be thresholded 
immediately. The output FM is placed in an array of row buffers of height 
N/2 that transform the rows and filter immediately the M£?'s produced. 
The MF output is placed in another vertical buffer where the traditional 
one-dimensional transform procedure is enacted. The FF output is placed 
in an array of row buffers identical to the original configuration, except that 
it is only N/2 tall. The same procedure is enacted on this data set that is 
half as small. Now this proceeds across the matrix in a similar manner as 
the one-dimensional case, except that the vertical buffers can be completely 
purged as the next column is reached. To count the number of elements that 
are necessary we can ignore the vertical buffers, which are subdominant. 
At the first level we note that there are N * (2L — 2) elements in the row 
buffer and (JV/2) * (2L - 2) in the FF output and FM output. Hereafter 
we will drop the 2L — 2 to simplify the counting since we are just looking for 
order. So we have N and N/2 and N/2. The FM output will proceed like 
the one-dimensional case. Therefore it will produce (N/2) log(iV) elements. 
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The FF out put will produce N/AFF and N/AFM. So we can see that 

the total number of elements will be N{\ + log(iV)) X)/=o ¥ ■ ^ ne sum * s a 
simple geometric sum that in the limit that J goes to infinity is bounded 
by 2. So the final tally of the necessary elements is 0(N\og(N)). The 
generalization to D dimensions is straightforward. One begins with a data 
structure of dimensions (N D ~ 1 )(2L — 2). One then performs a transform to 
produce two N D ~ 2 data structures. One performs a transform on these two 
structures to produce four N D ~ 3 structures. This process continues until the 
final transform where one has a single dimension. This transform is enacted 
and the M D elements are filtered. Appropriate lower dimensional transforms 
are applied to the mixed output, MMMFF ■ ■ ■ MF, FFFFFF ■ ■ ■ M, etc. 
The process is repeated for the FFF ■ ■ ■ FFF data set. In higher dimensions 
the algorithm becomes more complicated, but the idea is the same. And the 
leading order number of elements that need to be saved is O^N 0-1 log(iV)). 



10 Appendix I: Continuous Wavelets 

We begin by considering the continuous wavelet transform. The continuous 
wavelet transform is an alternate representation of a function, like a Fourier 
transform. Both continuous and discrete wavelets are built from a single 
function called a mother function. The notation, ip(x), is used to denote 
the mother function of a wavelet. 

Wavelets are built from translations and scale transformations of the 
mother function. Translations and scale transformations of ip(x) are defined 
by: 

iM*) := \ 8 \-*1>(—). (238) 
s 

The factor p is a parameter. The functions ipt,s(x) are the wavelets associ- 
ated with the mother function i[)(x). The wavelet ipt,s{x) has two continuous 
parameters. We investigate conditions on the mother function that allow one 
to expand any function in terms of wavelets. 

To choose the parameter p note that 



\s\- p ^( 



x — t. 



'I 



dx 



/oo 
\i){u)\ q du. (239) 
-oo 
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It follows that if p = 1/q the L^-norm of ip 



a: 



\ 1/9 

\^(u)\ q du) (240) 



is preserved under scale transformations. Thus for p = 1/q: 

IMI, = lhMI« for all s,t. (241) 

The continuous wavelet transform of / is defined by taking the scalar 
product of / with the wavelet ip s /. 

PCX) 

f(s, t) := / tf t {x)f{x)dx = (i/>,,u f) (242) 



s 
oo 



where asterik V indicates the complex conjugate for a complex mother func- 
tion. In what follows a / is used to indicate the wavelet transform of the 
function /. 

Parseval's identity for the Fourier transform implies that the wavelet 
transform can be expressed in terms of the original function and the mother 
function or alternatively in terms of their Fourier transforms: 

/(M) = (ik,*, /) = (&,*, /) (243) 
where the ~ indicates the Fourier transform defined by: 

^ 9 , t (k) = -= / e- ikx ^ 9 , t (x)dx (244) 



f(k) = I e- ikx f(x)dx. (245) 



2tt 

Note that Parseval's identity states (/, /) = (/, /). Using this with / = 
g + h and f = g + ih gives 

(g, g) + (h, h) + (g, h) + (h, g) = (g, g) + (h, h) + (g, h) + (h, g) (246) 

and 

(g, g) + (h, h) + i(g, h) - i(h, g) = (g, g) + (h, h) + i(g, h) - i(h, ~g) (247) 
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which, using the identities (g, g) = (g,g) and (h,h) = (h,h), gives the solu- 
tion to (gUD and (HZD: 

(g,h) = (g,h) (248) 

which is the form of Parseval's identity used in (J243j) . 

The Fourier transform of the wavelet ip s ^{x) can be expressed in terms of 
the Fourier transform of the mother function: 

1 r°° r — f 
i/> a>t (k) := -f= I e- lkx \s\- p ^j{- -)dx = 

-oo S 
■iksu —ikt I I — p+1 



/2tt 

1 



e~ lksu e- lkt \s\- p+1 ^(u)du 



oo 



| s |i-p e -<*ty( s fc). (249) 
The inner product of the Fourier transforms gives 

f(s,t) = $ a ,tj) = 



oo 



r s , t (k)f(k)dk 

) 

s\ l - p e lkt ij*(sk)f(k)dk. (250) 



Multiplying both sides of ()250|) by e %k ' 1 and integrating over t gives 



2tt 



' ' e- lk,t (^ t J)dt 



Is^risk^fik'), (251) 
where the representation of the delta function: 

1 f°° 

_ / e -i{ k '-k)t dt = 5{k , _ k y ( 252 ) 

was used to get (|251|) . 

The right-hand side of ()251|) is a product of the Fourier transform of the 
original function with another function. We can't divide by the function 
ip*(sk') because it might be zero for some values of k! . Instead, the trick is 
to eliminate it using the variable s. 
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Multiply both sides of this equation by ip(sk') and a yet to be determined 
weight function w(s) and integrate over s. This gives 

-| poo poo 

— w(s)ds dte- ik,t 4>(sk')f(s,t) = 

/oo 
w(s)rfs|s| 1 -V(sA; , )^(sA; , ) = f(k')Y(k') (253) 

where 

/•oo 

Y(k') = / rfsw(s)|s| 1 - p |^(sA; , )r. (254) 



In order to be able to extract the Fourier transform of the original func- 
tion, it is sufficient that Y(k') satisfies < A < Y(k') < B < oo for some 
numbers A and B. In this case 

i poo poo 

/(*) = 2^r(}b) i> w ^ ds J dte- M ij(sk)f(s,t). (255) 
It is convenient to rewrite this in terms of the wavelet basis: 

i poo poo 

~ f{k) = 2*Y(k)j w ^\ s \ P ~ lds J dti> s , t (k)f(s,t). (256) 
We define the dual wavelet by 

The dual wavelet is distinguished from the ordinary wavelet by having the 
parameters s, t appearing as superscripts rather than subscripts. 

The inversion formula can be expressed in terms of the dual wavelet by 

/oo poo 
w{s)\s\ p - l ds I dt^ s >\k)f(s,t). (258) 
OO J —00 

In order to recover the original function, take the inverse Fourier trans- 
form of this expressions: 



1 f°° 

f(x) = -j=J dke lkx f(k) 
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oo 

p-1. 



oc 



w(s)|s| p -Ms / dt^ t (x)f(s,t) (259) 



where 



^\x) = -= \ dke lkx ^\k). (260) 



OG 



In general this is a tedious procedure because the dual wavelet ip s,t (x) must 
be computed using (j257J) and ()260|) for each value of s and t. If the dual 
wavelet also had a mother function, then it would only be necessary to Fourier 
transform the "dual mother" and then all of the other Fourier transforms 
could be expressed in terms of the transform of the "dual mother" . 

The first step in constructing a "dual mother" is to investigate the struc- 
ture of the dual wavelets in x-space: 



iP s ' f (x) = -= \ dke lkx ^ s \k) 

V27T J -„o 



1 



where 



_ , dke lkx ——\s\ 1 - p e- lkt ij(sk) 
^ s '°(x - 1) 



If 00 l 
r>\x) = -= / dke^——\s\ 1 ~^(sk). 



This shows for a single scale the dual wavelet and its translation can be 
expressed in terms of a single function. This is not necessarily true for the 
dual wavelet and the scaled quantity. 

1 f°° 1 
= -= / dke^——\s\^(sk) 
V27ri-oo 27rr(fc) 



1 



due iu « — \s\- p ^{u). 



^ttJ-^ 2ttY(u/s) 

This fails to be a rescaling of a single function because of the s dependence 
in the quantity Y(u). It follows that if a weight function w(s) is chosen so 
Y(u/s) = Y is constant, the dual wavelet will satisfy 

1 f°° 1 

ij; s ' {x) = -= due iu s——\ s \-P^{u) = \s\- p ^{x/s). (261) 



2W-oo 27TF 1 
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In this case Y(u) is a constant which we denote by Y. The function ip 1,0 (x) 
serves as the dual mother wavelet. 
To determine w(s) note that 



Y(sk) = / dtw(t)\t\^ p \^(tsk)\' 

Let t' = st to get 

/oo 
dtiu(*)l*M^(tsfc)| 2 
-oo 

/oo 
rft'w(t7s)|t| ,1 - p |^(t , A;)p 
-oo 



IP-2 

./ — oo 

This will equal Y(k) provided 



w {t') = \s\ p - 2 w(t'/s) or w(s) = \s\ p - 2 w(l). 
With this choice 



° (II - 2 



= y = «,(i) / -|^)| 

./-oo l r l 

Assuming this choice of weight the admissibility condition becomes 

< A<Y < B < oo. 

Having computed the constant Y it is now possible to write down an 
explicit expression for the dual mother wavelet: 

if, a >°(x-t) = 

\s\- p r , ■■■■(« " 1 



Letting k = u/s 



L due ' sr* u) 



/OO -1 



27T ^ - v 

/•oo 1 

ifca;. 
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dk—e lkx ^ t {k). 



This has the form 

r>\x) = ^ Mx) - (262) 

Thus the inversion procedure can be summarized by the formulas: 

/oo pco 
\s\ 2p - 3 ds / dtij s '\x)f(s,t) (263) 
-oo J —oo 

r°° dt ~ 

Y= -Mt)\ 2 (264) 

J-oo \ l \ 



r\x) = (265) 

i> s ,t = \s\-^{^). (266) 

The mother function must satisfy < Y < oo. This requires that the 
Fourier transform of the mother function vanishes at the origin. This is 
equivalent to saying that the integral of the mother function is zero. 

Using the representation for the wavelet transform gives a representation 
of a delta function: 



5(x-y) = 

/oo poo 
\ S \ 2p - 3 ds / dt^\xW s M = 
-oo J — oo 

> poo 

\ s \ 2 ^ 3 ds / dt^, t (x)r s M- 



2nY 



We can also use this representation of the delta function to formulate a 
Parseval's identity for continuous wavelets 

i poo poo 

VJ) = ^yJ Js\ 2p ~ 3 dsj Jt\f(s,t)\ 2 . (267) 
Consider the example of the Mexican hat wavelet. The mother function 

is 

ij(x) = -^{x 2 - l)e- x2/2 . 
v 2tt 
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To work with the Mexican hat mother function it is useful to use prop- 
erties of Gaussian integrals: 



oo 

e -ax*+bx+c dx 



oo 



Change variables to y — y/a(x — j-) to obtain: 

2 

e~ y dy = 



2a 

b 2 



e 4a . 



This can be used to compute the Fourier transform of the Mexican hat 
mother function: 



i r°° 

^(fc) = -= / e- ikx i)(x)dx = 
V2tc J-oo 



I r°° 

II '2 



2tt /_ 



x~ - l )e- x l 2 ~ lkx dx. 



-oo 

To do the integral insert a parameter a which will be set to 1 at the end of 
the calculation: 



(_ 2 A_1)J_ r e -x*a/2-ik Xdx = 
da '27T./-00 

. n d „ . 1 /2tt fc 2 
-2— - — \ — e - ^ = 



( 2~ l )\h^ e 2a - 

a er V 2na 
In the limit that a — > 1 this becomes 

— k 2 e~^~. 
2tt 

Using this expression it is possible to calculate the coefficient Y 




/°° dh ~ 
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1 

2^ 



\k\ 3 dke- k2 



-I 

k Jo 



oo 
oo 



k 3 dke~ k2 



Inserting a parameter a which will eventually be set to 1 gives 

1 d C°° 
Y = — (-- f ) / 2kdke~ ak2 = 
2vr da J J 

2ir da a J 
1 

2tt' 

This satisfies the essential inequality < Y < oo which ensures the admissi- 
bility of the Mexican hat mother function. 

The expression for the wavelet transform and its inverse can be written 

as: 

fat) = \sr r ) 2 - l)e-^^f(x) = 

J-oo y/Zrr s 

/OO 1 
du^={u 2 -l)e- u2/2 f(su + t). 
-oo V 27T 

where x = su + t 

The inverse is formally given by 

/(*) = Jjs\ 2p - 3 ds J^dt^f(8,t) = 

\s\ 2p - 3 ds r <ft_L| a |-P((^*)' - l)e-^) 2 /V>,0 = 
J-oo J-oo V 2ir s 

i = r \*r*d8 r dt{{ x -^f - i)e-^f(s,t) = 

2% J-oo J-oo S 

-i i>oo poo 

—= \s\ p - 3 ds du(u 2 - l)e^ 2 / 2 /(s, SU + x) 

V27T J-oo J-oo 
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where t = su + x. 

Initially we were concerned because we were representing an arbitrary 
function by a linear superposition of functions that all have zero integral. 
We could not understand how wavelets could be used to represent a function 
with non-zero integral. 

We tested this by computing the wavelet transform and its inverse for 
a Gaussian function with the Mexican hat wavelet. The original Gaussian 
function was recovered. 

The resolution of this paradox has to do with the difference between L 1 
and L 2 convergence. The wavelet transform has a vanishing L 1 norm, but 
the L 2 norm is non-zero. 



11 Appendix II - Spline Wavelets 

We use the convention which defines the Fourier transform of a function f(x) 
as 

/oo 
e~ ikx f(x) dx , (268) 
■oo 

and the inverse transform by 

POO 

2tt 



f( x ) = — I F(k)e ikx dk. (269) 



For this convention, Parseval relation is 

/OO -I /'OO 

f*(x)g(x) dx = — / F*(k)G(k)dk. (270) 
■oo 27T > /_ 00 

The cardinal B-splines, N m (x), are defined by first defining 

( 0, if x < 
N^x) = < 1, if < x < 1 (271) 
1 0, otherwise. 

Then for m > 2, N m (x) is defined recursively by the convolution integral 

/oo 
N^x-^N^dt 
-oo 

1 

Nr^x-^dt. (272) 
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Since N m (x) is denned by a convolution integral, the Fourier transform will 
be defined by a product. To show this, we evaluate the Fourier transform 



~ lkx N m (x) dx 

"OO 



Nm(k) = [ 

J —oo 

/oo 
dxe~ ikx I Nn-^x-ijN^ijdt 
-oo J —oo 

/OO POO 
dte-^N^t) I e- lk{x - t] N m ^(x-t)dx. (273) 
-oo J —oo 

Now setting x — t = y in the second integral, we find 

/oo poo 
dte-^N^t) / e-^N^Mdy 
-oo J —oo 

POO 

= N m -i(k) / e'^N^dt 



- _ im 

= Ni(fc) , 



(274) 



where 



iVi(A;) 

■oo 

1 



-I 



r lhx dx 



= e 



ik/2 S H k / 2 ) 

k/2 



(275) 



For this example we use the quadratic spline shifted to the left by one 



unit 



B(k) = e lk N 3 (k) 



-ik/2 



"sin(fc/2)" 
k/2 



(276) 
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Now evaluating the inverse transform using Maple, we get 

fO, ifx<-l 
\(x + lf, if-l<x<0 



B(x) 



\-{x-\) 2 , ifO<x<l 



(277) 



|(^-2) 2 , 
10, 



if 1 < x < 2 
otherwise. 



The splines are not orthogonal; however, we can use them to construct a 
scaling function <j)(x) which has the ort honor mality property 



/oo 
(f)*(x — l)(f)(x — m) dx = 5i m ■ 
-oo 



(278) 



To do this we follow the procedure given in the books by Chui 1 and Daubechies 2 . 
Note that this a general procedure; we are using the spline as a convenient 
example. The method gives an expression for the Fourier transform, (j>(k), of 
4>(x). The Fourier transform of <f>(x — I) is given by 



(*) = f 

j — t 



e~ lkx <\){x - I) dx 



= e 



-ikl 



-ik(x—l) 



4>(x — I) dx 



^—ikl 



Now we show that if 



i r 



ikm 



4>(k) dk = S mfi , 



(279) 



(280) 



then the functions are orthogonal. To show this, we use the Parseval relation 



J — c 



<p*(x — l)(f)(x — m) dx = — 

Z7l 



2vr/ 



■oo 



(k)4> m {k) dk 
l 4>*{k)e~ ikm 4){k) dk 



1 Charles K. Chui, An Introduction to Wavelets, Academic Press, 1992 
2 Ingrid Daubechies, Ten Lectures on Wavelets, SIAM, 1992 
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1 

2^ 

Oj-m,0 
film ■ 



m) 



dk 



Finally, we show that if we can find a 4>(k) such that 



\(p(k + 2nl) 

l=—oo 



(281) 



(282) 



then the functions are orthonormal. The infinite sum in Equation ()282j) is 
periodic in k with a period of 2tt; thus it has the Fourier series expansion 



\4>(k + 27il) ~ = c 3 eik3 

l=—oo j=—co 

where the expansion coefficients are give by 

»2tt oo i 2 



(283) 



_ J e -*h J- \^k + 2nl) 



-E 

2tt ^ 



oo r i-K 



l=— oo 



-* fci U(Jb + 27rZ) 



j 00 p2tt(1+V) 



E 

i=— 00 

/>00 



2tt 



2tt 



27ii 

(Ay 



-i(k-2wl)j 



dk . 



dk 



(284) 



Since the sum in Equation ()282|) is equal to one, Cj = fij t0 , and one finds 



1 

2^ 



-ikj 



4>(k) dk = fijfi . 



(285) 



Thus, the functions are orthonormal. 

Now given a function, B(x), we construct a scaling function by taking its 
Fourier transform and defining 



}(k) 



\B{k + 27il) 



(286) 



j=— 00 
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This function satisfies Equation (|282j) . and the <fi{x) will have the orthonor- 
mality property given in Equation (J278|) . To evaluate the infinite sum in 
Equation ()286j) . we use the finite Fourier series expansion of the function 



g(k)= J2 \B(k + 2nl) 



l=—oo 



(287) 



This function has period 2ir, and the Fourier expansion has the form 



Cj e 



ijk 



j=-oo 



Following the derivation in Equation 
given by 



(288) 

the expansion coefficients are 



2- 



— J e^ k g(k)dk 
1 r 2lT 00 

J U l=-oo 

— / e- ijk \B{k)\ dk 
2vr J_ 00 I I 

— / B*(k)e~ ijk B(k)dk 

/oo 
B*(x)B(x - j)dx, 
■oo 



(289) 



where the Parseval relation was used for the last step. The integral in Equa- 
tion ()289|) is easy to evaluate for the B-splines, and we find 



f -i- if ? 

120 ' •> 

60 ' 
11 
20 

M> if -? = 1 

T20' if i = 2 



B*(x)B(x - j)dx = < 



, 0, otherwise. 



(290) 
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Using these coefficients for the expansion given in Equation (|288|) . we find 

(291) 



9 ^ = 27) + 1! C ° S ^ + 60 cos( ^ 2A; ) • 



To find 4>(x) the inverse Fourier transform of <f)(k) must be done numer- 
ically; however, there is a nice method which gives an efficient algorithm. 
Since the function g{k) has period 2tt, we can use the expansion 



1 



E 



-ink 



where the coefficients 



2tt 



7t ginfc 



dk 



(292) 



(293) 



must be computed numerically. For the B-spline, the g(k) is an even function 
of k, and one finds 

1 ^ C °< nk) dk. (294) 



In addition, from Equation (|294J) we see that c_ n = c n . Now, using Equation 
(jH2D we get 



0(z) 



1 

2^ 
1 

2^ 



— oo 
oo 



B(k) 



-ink 



l 

2^ 



oo 

E 

— c 



n=— oo 
oo 



E 

n=— oo 
oo 

= c n 5(x - n) . 

n=— oo 

Now we need to find the wavelet ip(x) with the properties 

ip*(x — l)cj){x — m) dx = , 



(295) 



(296) 
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and 

/oo 
if>*(x-l)ij)(x-m)dx = 6i m . (297) 
-oo 

To do this we introduce the functions 

(f>-i,i(x) = V2<J)(2x - I) . (298) 
The Fourier transform of these functions is given by 

/oo 
e- ikx <p-^i{x)dx 
-oo 

= V2 e~ ikx (j)(2x - 1) dx (299) 



■oo 



and setting 2x — I = y yields 



4>-iM = 77f~" /2 I e- tky/2 Hy)dy 



= —e- ikl ' 2 4>{k/2). (300) 
v2 

Using the 0_i jra (x) as an orthonormal basis set, we can write 

oo 

= E h n<l>-l,n{ X ) ' ( 301 ) 
n=— oo 

where 



<j>U„ n {x)<t>{x)dx.. (302) 

This is the scaling equation for this system. In this case there are an infinite 
number of non-zero scaling coefficients. Since the <fi(x — n) are orthonormal, 
the h n must have the property 

IM 2 = 1- (303) 

n=— oo 

Taking the Fourier transform of Equation (j301|) gives 

1 oo 

= ^ E h n e-^ 2 4>(k/2) , (304) 
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which can be written as 



4>{k)=m {k/2)4>{k/2), 



where 



-ikn 



(305) 
(306) 



Using Equations ()282|) and (j305|) we see that 



Y \4>(2k + 2nl) 



l=— oo 



2 _ 



|mo(A; + 7r/) <j){k + -nl) 

=— oo 

El I 2 I ~ 

mo(fc + 7r/) </>(&; + 7r/) 

,ei>en 

+ Y \ m o{k + ?rZ) </>(£; + ttZ) 

°° 2 

m (A; + 27r/) 4>(k + 2nl) 



2 _ 



l=— oo 



2 _ 



+ |m (A; + 27r/ + 7r) 0(fc + 2vr/ + vr) 

l=— oo 
00 1 2 

|m (fc)| 2 ^ U(A; + 2ttO 



!=— oo 



+ \m (k + vr)| 2 ^ + vr + 2nl) 



l=— oo 



|m (/c)| 2 + |m (fc + 7r) J" 



(307) 



where we have used the periodicity of m (k). The sum on the left-hand side 
of Equation (j307J) is equal to unity; thus, we have shown that 



\m (k)\ 2 + \m (k + n)[ 



(308) 



Now we use a similar procedure to find ip(x). Using the 0_i jn (x) as an 
orthonormal basis set, we write 



x 



(309) 
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with 

/oo 
<t>*_ 1)n (x)1>(x)dx. (310) 
-oo 

Taking the Fourier transform of Equation (|3U9|) gives 

1 oo 

^) = ~m E /ne-^ 2 0(A;/2) 

v n=— oo 

= m 1 (k/2)4>(k/2), (311) 

where 



oc 



v n=— oo 

If mi(k) has the same property as that given in Equation (|3()8|) for mo(k), 
then the functions ^(x — m) will be orthonormal. In addition, we want 
ip(x — n) to be orthogonal to <f)(x — m). Thus, we want to find a mi(k) such 
that 



30 1 /"OO 

^*{x-n)(j){x-m) = — e ikn ^*(k)e~ ikm ^(k)dk 

oo ^7T J-oo 



OG 



— / e i{n - m)k ^*(k)0(k)dk 

2tt 



- / dke i(n - m)k ^*( k + 2vr/)0(fc + 2ttZ) 

^° l=-oo 



2tt 

= 0. (313) 
This condition is satisfied if 

oo 

^2 + 2irl)4>{k + 2-kI) = . (314) 

(=— oo 

Substituting Equations (|3U5|) and (|311|) into Equation (|314|) and replacing k 
by 2k gives 

oo 

Y m*i{k)4>*{k + 7rl)m {k)4){k + ttI) = . (315) 

!=— oo 

Regrouping the sums for odd and even I, and following the procedure used 
in Equation ()307|) gives 

m[(k)m Q (k) + m\{k + 7r)m (fc + tt) = . (316) 
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This condition will be satisfied if we choose 

mi(k) = e- ik m* (k + vr) . (317) 

Note this choice for mi(k) is not unique; we can multiply m\{k) by any func- 
tion p{k) which has period n and \p(k)\ = 1, and still satisfy the constraints 
on mi{k). Substituting this result into Equation (|311|) gives 

^(k) = e- ik/2 m*{k/2 + 7c)4){k/2) 

! 1 E Ke l{k/2 ^ )n kk/2) 



e -ik/2 



V2. 

oo 

-7E E {-l) n Ke-^ n+1)k/2 Kk/2) 



n=— oo 

oo 



n=— oo 
oo 



V2 : 

oo 

= E (-l)X^iW ■ (318) 

n=— oo 

Replacing n by — n + 1 in Equation ()318|) gives 

oo 

m = - e (- i ) B ^i?-i,n(*) • (319) 

n=— oo 

For convenience, we drop the minus sign in front of the sum, and write 

oo 

$(k) = E 9n4>-i,n{k), (320) 

n=— oo 

where g n = (— l) n h_ n+ i for /i n a real number. Taking the Fourier transform 
of Equation (|32()jl gives the result 

oo 

^( X ) = E 9n<l>-l,n(x) 
n=— oo 

oo 

= >/2 E ^^( 2a; - n ) ■ (321) 

n=— oo 

To evaluate ip(x) we use the expansion given in Equation (|295jl in Equa- 
tion (j321|) . This gives 

oo oo 

if)(x) = V2 ^2 g n E c mB{2x -n-m). (322) 



n=— oo m=— oo 
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Now replace m by I — n, this gives 



^2 ^ 9nCl- n 



l=—oo 

00 

= d iB(2x - /) . 

l=—oo 

Numerical Methods 



5(22; - I) 



(323) 



To determine the h n we need to evaluate the overlap integral given in 
Equation (13021 . From Equations (EMI) and we find 



-i,n(aO = c m B(2x - m-n) 



m=—oo 



Then using 



S(x)= £ b j B{2x-j) 



j=-oo 



where, the 6j for the quadratic B-splines are given by 



(324) 



(325) 



we can write 



/ 1 

4' 


if j = -l 


3 
4' 


if 3 = 


< 2 

4' 


if j = l 


1 

4' 


if j = 2 


,0, 


otherwise. 


3 


00 



£ c n £ b 3 B(2x-2n-j) 



n=—oo j=—oc 



2J 2^ C «^~2n 

l=—oo [n=—oo 
00 

£ s,B(2a; - Z) . 
/=— 00 



5(2a; - /) 



(326) 



(327) 
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Using these expansions, the overlap integral is given by 



OO OO 

V2 CmSl B{2x-m-n)B{2x-l)dx 

m=-ool=-oo J-oo 
OO OO ^ /»oo 

\P1 c m s\- j B(x — m — n + l)B(x) dx 

m=-oo Z=-oo J-oo 
.j oo 2 /»oo 

— ^ ^ Y CmSm+n^j j B{x)B{x - j) dx , (328) 



m=— oo j=—2 



where, we have set / = m + n — j in the second summation. The values for 
the integrals of the quadratic B-splines are given in Equation (|289p . 

To derive Equation ()325|) . we use sin(26 l ) = 2cos(#) sin(#) to write Equa- 
tion (j2Zn)) as 



B(k) 



-ik/2 



-ik/4 



-ik/A 



cos(/c/4) sin(fc/4) 
kjl 

e ifc/4 _|_ e -ik/4 



D ik/4 _|_ g— ik/4 



-ik/2 



sin(fc/4) 
fc/4 



B{k/2) 



e ik/2 + 3 + 3 e -<*/2 + e -« 



8 



B{k/2). 



(329) 



Now taking the Fourier transform of (|329|) and using 

/oo -t poo 

e~ ikx B{2x- j)dx = - e ik]/2 e~ ikx/2 B{x) dx 
-oo ^ J —oo 



-e ikl,2 B(k/2) 



(330) 



we find 



5(x) = ^B{2x - 1) + ^5(2a;) + ^B{2x + 1) + ^5(2a; - 2) . (331) 



The authors would like to thank Mr. Fatih Bulut for a careful proof 
reading of these notes. 
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